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Abstract 



I consider the initial-boundary-value-problem of nonlinear general relativistic vacuum space- 
times, which today cannot yet be evolved numerically in a satisfactory manner. Specifically, I 
look at gauge conditions, classifying them into gauge evolution conditions and gauge fixing con- 
ditions. In this terminology, a gauge fixing condition is a condition that removes all gauge degrees 
of freedom from a system, whereas a gauge evolution condition determines only the time evolu- 
tion of the gauge condition, while the gauge condition itself remains unspecified. 1 find that most 
of today's gauge conditions are only gauge evolution conditions. 

1 present a system of evolution equations containing a gauge fixing condition, and describe 
an efficient numerical implementation using constrained evolution. 1 examine the numerical be- 
haviour of this system for several test problems, such as linear gravitational waves or nonlinear 
gauge waves. 1 find that the system is robustly stable and second-order convergent. 1 then apply it 
to more realistic configurations, such as Brill waves or single black holes, where the system is also 
stable and accurate. 
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1 Introduction 



I want to simulate spacetimes containing black holes. 

1 am not alone with this project; it has a long history. This was the goal of the Binary Black 
Hole Coalescence Grand Challenge in the USA, and is also part of the programme of the SFB 382 
in Germany which pays my salary. One wants to numerically simulate spacetimes that contain 
singularities, where the spacetimes are time-dependent and have no symmetries that reduce the 
dimensionality, and where no approximations are possible to simplify the nonlinear equations. 

Such calculations are often justified by referring to gravitational wave detectors such as GEO 
600, LIGO, or LISA, which will soon produce data which then need to be compared to theoretical 
calculations. 1, however, think that the ability to and the knowledge how to simulate fully generic 
spacetimes have value in itself and need no further justification. 

The two major problems that one encounters today when simulating spacetimes are the com- 
putational expense and numerical instabilities. Three-dimensional time-dependent simulations 
always produce much data, and always need supercomputers to run on. Supercomputers are 
awkward to use, and it is time-consuming to experiment and try new ideas. 

The numerical instabilities that are encountered in spacetime simulations are a problem of a 
different kind. For many years now, simulations have been plagued by instabilities of unknown 
origin. Simulations run fine for some time, but in many cases they encounter an unphysical and 
poorly understood unboimded growth. 

Many approaches have been tried to find a cure for these instabilities. People have blamed the 
formulation of the equations, the boundary conditions, the gauge conditions, and the discreti- 
sation methods. It is clear that errors in any of these can cause a simulation to fail. The only 
successful methods to remove the instabilities known today rely on a combination of several steps 
that seem to be random, looking like a magician's recipe, and can only be justified after the fact. 
While it is now possible e.g. to evolve a static black hole for a long time, or to run a black hole 
collision for a reasonable time, the overall situation is still far from satisfactory (see e.g. |SY02|). 

In this thesis 1 look at gauge conditions, which are one of the possible causes of instabilities 
mentioned above. 1 distinguish between gauge fixing and gauge evolution conditions. Gauge fixing 
conditions are "real" gauge conditions, whereas the weaker gauge evolution conditions describe 
only the time evolution of an (unspecified) gauge condition. The majority of today's gauge con- 
ditions are only gauge evolution conditions in this terminology. I propose a set of gauge fixing 
conditions, and show that their additional power can lead to a stable evolution, even without 
using any magic ingredients. 

While experimenting with different formalisms, 1 am willing to sacrifice speed. 1 do believe 
that it is enough to have a formulation that can at least in principle be implemented efficiently. 
Computers become faster every year, and run time on supercomputers is rather easy to obtain 
for research purposes. It is better to start with a slow, but stable formulation, and optimise later 
on, than trying to find stability among fast implementations. It is also better to have a simple 
scheme for doing things, than having to follow a long recipe containing many steps, each of which 
is necessary. 
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1 Introduction 

This thesis is structured as follows. I start with a general introduction to simulating vacuum 
spacetimes and the problems encountered therein (chapter |3. I then explain gauge fixing condi- 
tions in general and the condition I chose in particular (chapter^, as well as a system of evolution 
equations that is suited to this kind of condition (chapter |4j. After that I discuss some bound- 
ary conditions (chapter |5j and initial data (chapter |6j that I use later on. I briefly mention some 
facts about the code that I developed (chapter O and move on to present in some detail the test 
cases and their results that demonstrate the properties of my gauge fixing conditions (chapter 
I make some concluding remarks (chapter|9)/ and add some appendices containing equations and 
explanations that are likely not of general interest. 
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2 Simulating spacetimes 



A spacetime can be characterised by its four-metric g,j,v, which has to satisfy the Einstein equation 
G^v = SnTj^^. I will restrict myself in this thesis to the vacuum case where 7^^-^ = 0. (The notational 
conventions are explained in appendix ICSI ) The four-formalism couples the spatial and temporal 
degrees of freedom. The astrophysical interpretation, however, calls for an initial-boundary-value 
formulation. 

2.1 The ADM formalism 

One commonly used way of transforming the Einstein equations into such a system of time evo- 
lution equations is the so-called ADM formalism IADM62I . This is one form of a so-called 3-1-1 
split into 3 spatial and 1 temporal dimensions. It leads to the primary variables three-metric y,y 
and extrinsic curvature K^j. The time evolution equations contain the freely specifiable quantities 
lapse a and shift 13'. The three-metric and extrinsic curvature have to satisfy two constraints, the 
Hamiltonian constraint H and the momentum constraint M, . 

In order to evolve a spacetime in an ADM-like formalism, one needs initial data, boundary 
conditions, and a choice of lapse and shift. Additionally, one can choose to enforce the constraints 
or a gauge. For reference, the ADM variables are given in appendix lA. 1.11 the ADM evolution 
equations in appendix lA. 1.21 and the ADM constraints in appendix lA.1.41 

At first sight, the primary variables of an ADM system (three-metric yjj and extrinsic curvature 
Kij) seem to contain six physical degrees of freedom, i.e. twelve independent quantities. However, 
only four of those twelve have a true physical meaning. Out of the twelve apparent independent 
quantities, the constraints eliminate four, and the choice of gauge eliminates another four. Thus 
there are only four independent quantities left, corresponding to the two degrees of freedom for 
gravitational waves. 

These days^, the original ADM system has gone out of fashion in the numerical relativity com- 
munity. Instead, one uses extended ADM-like systems, often called ctADM ("conformal traceless 
ADM") systems. They have two more variables, a conformal factor xp for the three-metric and the 
trace K of the extrinsic curvature. The three-metric is replaced by a conformal three-metric y,y, and 
the extrinsic curvature by a conformally rescaled traceless extrinsic curvature Ajj. These systems 
also have two additional constraints, namely det y,j = 1 and trace Ajj = 0. 

Both the ADM and the ctADM systems exist in multiple variations. The constraint equations 
are of the form C = 0, and it is therefore possible to add arbitrary multiples of the constraints C 
to the right hand sides of the time evolution equations. Doing this can dramatically change the 
behaviour of a system. 

The ctADM systems have the advantages that the new variables xp and K have a direct inter- 
pretation as constraint and gauge quantities, ip is closely related to the Hamiltonian constraint H, 

write this in 2002 
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2 Simulating spacetimes 

while K is related to the lapse a. These relations make it easier to influence the time evolution 
of these quantities. It also makes it potentially easier to enforce the constraints with the York- 
Lichnerowicz method, which also relies on a conformal traceless decomposition. 

Both the M aya |Pen LSOH lKLL+oTl and the AEl variants of the BSSN systems INOK87IIN0891 
ISN951 IBS99I are variants of ctADM systems, where certain derivatives of the conformal three- 
metric have additionally been made independent variables, which leads to three additional con- 
straints. The additional variables capture gauge degrees of freedom that are closely related to the 
shift, and thus allow further control over the system. 

2.2 Constrained evolution 

Sometimes it is necessary or convenient to enforce the constraints explicitly. This can either be a 
part of a scheme to obtain initial data, or it can be part of an evolution scheme, e.g. to prevent 
deviation from the constraint hypersurface due to accumulated numerical errors. Enforcing the 
constraints can also help stability, as e.g. described by Evans |.Eva89 1 or Marsha and Choptuik 
IMC96I . The constraint equations contain no time derivatives, and also do not contain the lapse or 
the shift. They connect the three-metric and extrinsic curvature only. 

When spherical or axial symmetry is assumed, then some constraint equations are algebraic 
conditions for some of the evolved variables. That means that these constraints can be enforced 
rather easily. This is unfortunately not the case in three dimensions. There, the constraint equa- 
tions are not suitable for enforcing in their ADM form. However, they can be rewritten so that they 
are elliptic equations in certain quantities that are connected to the three-metric and the extrinsic 
curvature. In the York-Lichnerowicz formalism IYor73l lYor74l ICoo02l section 2.2.1], a ctADM 
system is extended to include a vector potential for the conformal traceless extrinsic curvature. 
The constraints determine the conformal factor and this vector potential. This requires variable 
transformations before and after solving. 

As it is somewhat expensive to solve elliptic equations, it has become customary to perform 
unconstrained evolution for three-dimensional simulations, i.e. to evolve without regard to the con- 
straints, monitoring the constraints only as a measure of quality for the time evolution. Given ini- 
tial data that satisfy the constraints, the Bianchi identities then guarantee that the constraints will 
be satisfied at all times, of course only modulo numerical errors. Choptuik ICho91l argues that 
the order of accuracy of an unconstrained evolution will not be worse than that of a constrained 
evolution. 



2.3 Picturing spacetimes 

There are two fundamentally different points of view that one can take to interpret a spacetime. In 
the following, 1 use the terms grid points and coordinate system in an equivalent manner. One can 
think of coordinate lines forming a grid, and the intersections of these coordinate lines are then 
the grid points. 

The coordinate based view treats the primary variables similar to hydrodynamic or electrody- 
namic field quantities. One pictures a non-physical background grid, containing regularly 
spaced grid points forming the numerical grid. On each grid point, or in each grid cell, live 
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2.3 Picturing spacetimes 



the primary variables, i.e. the three-metric and extrinsic curvature. The background grid has 
nothing to do with the evolved spacetime, it exists only in a flat coordinate space. 

Time evolution consists of updating the evolved quantities through source and advection 
terms. The boundary of the grid stays in a fixed place in the coordinate space. This picture 
is the one that is most often used when quantities are plotted or displayed. Instabilities 
show up as certain metric coefficients becoming negative or some quantity growing without 
bound. In this picture, the apparent horizon of a static black hole can move across the grid. 

In the physically based view, one looks at the location of the coordinate grid points in the phys- 
ical spacetime. It is not always easy to picture (let alone draw or plot) a curved spacetime, 
but much of general relativity has to do only with strange coordinate systems, and applies 
even in a flat spacetime. It is instructive to use this picture to look at a time evolution of flat 
space in a nontrivial, time-dependent coordinate system. 

In this picture, time evolution consists of choosing the locations of the next layer of coordi- 
nate grid points. That is, it is the grid that is constructed layer by layer, while the (unknown) 
spacetime is not evolved, but only discovered. Instabilities show up as grid points getting 
too close or too far apart, grid cells becoming too elongated, or grid lines crossing. In this 
picture, the apparent horizon of a static black hole stays fixed, while the grid points can move 
across it. 

The coordinate based view is much easier to draw or display, and is hence what is usually 
presented. The physically based view is much more difficult to display, and requires in general 
some (artificial) embedding. However, the physically based point of view makes it also much 
easier to interpret gauge modes. For gauge modes, the embedding stays the same, only the grid 
points move. 

As an example, consider the time-dependent metric 

ds^ = -7-. df + 2 dtdx (2.1) 



+ 



2 



{2Atx) 

1 + 



dx^ + dy^ + dz^ 



with the free parameter A. For A = 0, this metric is identical to the Minkowski metric. This metric 
happens to describe a flat spacetime for all values of A, which is not obvious to see. Figure IZTl 
shows the graphs of the three-metric, extrinsic curvature, lapse, and shift for the case A = 1/2. 
This figure uses the coordinate based view 

Assume now that this deviation from the Minkowski metric is due to an unwanted gauge mode. 
It is possible to arrive at the conclusion that this metric is identical to the Minkowski metric after 
a gauge transformation by looking at the field quantities shown in fieure lZTl Yet it is not trivial to 
do so, nor can the nature of the gauge transformation (i.e. the gauge mode) easily be found out. 
However, one wants to know what gauge transformation was applied so that the gauge mode can 
be counteracted in a numerical evolution. 

Figure shows the same metric, but now uses the physically based view. That is, shown are 
the grid points of the coordinate system described by the metric l!2.1> , using a Cartesian coordinate 
system (i.e. the Minkowski metric) as reference. This graph is only possible because the spacetime 
is flat; otherwise one would need an embedding. It is obvious from this graph that the gauge 
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2 Simulating spacetimes 





Figure 2.1: Evolved field quantities vs. time and space 
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2.4 Interpreting the quantities 
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Figure 2.2: Location of the coordinate grid points in the (flat) physical spacetime 



transformation is a very simple one, namely just t' = f (1 + Ax'^) with A = 1/2. Knowing this, 
one could e.g. choose a different lapse to prevent the time slices from curling up. While this can 
in principle also be deduced by looking only at the primary field variables, it is certainly easier in 
the physically based view. 



2.4 Interpreting the quantities 

It is quite useful to have a mental picture of what the quantities in an ADM or ctADM system 
mean: 

• The conformal factor determines the volume that is contained in a grid cell; the volume is 
given by y^det y,y = xp^. 

• The diagonal elements of the three-metric determine the physical length of the coordinate 
lines between grid points within a time sUce. The off -diagonal elements determine the angles 
between the coordinate lines. 

• The trace of the extrinsic curvature determines the mean curvature of the layer of grid points 
forming the current time slice. It is closely related to the time derivative of the conformal 
factor. It is also closely connected to the lapse, as it is mostly the lapse that determines the 
change of the mean curvature in time. 

• I have foimd no immediately useful interpretation of the extrinsic curvature tensor. One can 
view it as the time derivative of the three-metric. 

• The lapse is the physical distance between two time slices, measured in proper time of an 
observer that is at rest in those slices. It is thus the ratio between proper and coordinate time 
for such an observer. When the lapse is zero, coordinate time increases, but no proper time 
elapses. 

• The shift determines the difference between the spatial coordinate systems of two neighbour- 
ing time slices. It is the velocity of the coordinate system with regard to a resting observer. 
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2 Simulating spacetimes 

Given the three-metric and extrinsic curvature of a time slice, there is no direct way of knowing 
just where in a given spacetime the current time slice is located. One has to work out how the time 
slice is embedded in the spacetime, and then compare the resulting three-metrics and extrinsic 
curvatures. This is in general a difficult matching process, and 1 know of no attempts to try this 
in four dimensions (but see the Lazarus project IBCL02I '). An efficient algorithm for this problem 
would surely be very useful for comparing results of runs with different formulations, or of runs 
with different coordinate conditions LSal021 lHaw02l |3pp| ■ 

2.5 Instabilities 

Three-dimensional numerical simulations of nonlinear general relativistic spacetimes have been 
plagued by instabilities. To date, it is still safe to say that these instabilities are not well understood. 
There are many ideas about on how to remove these instabilities, but none have truly succeeded 
so far. 

I use the term instability here in a rather loose manner. I take it to mean that a certain quantity 
grows exponentially (or faster) and without bound. 1 want to use this term here not for a theoretical 
argument, but only to describe empirical results from simulation runs. 

1 do not want to use the definition of instability from the theory of hyperbolic systems. 1 think 
that this definition is not too useful in this case, as it considers an exponential growth to be stable. 
(On the other hand, it introduces the concept of resolution independence, which is crucial for 
numerical simulations.) 

Jeffrey Winicour has introduced the very practical definition of robust stability ISGB W00IISSW02I . 
which allows for at most polynomial growth even when arbitrary noise is continuously injected 
into the system. This can easily be tested in a numerical implementation. 

Apart from being stable, a system must also be convergent. That is, it must be able to track 
a given solution arbitrarily closely, where the deviation should depend on the resolution e.g. in a 
quadratic manner. 1 will concentrate on stability, rather than convergence, in the following section. 
That does, of course, not mean that 1 believe convergence to be unimportant. 

Let me try to classify the possible instabilities of a system implementing the Einstein equations 
in a 3 + 1 split: 

Physical instabilities have their cause in the spacetime itself. A spacetime can contain singulari- 
ties, and if the part of the spacetime that is described by the numerical metric gets close to 
or contains such a singularity, there will be a physical instability. There is not much one can 
do against such a singularity. Possible approaches include choosing a slicing that does not 
intersect the singularity, or excising that part of the spacetime, or trying to "subtract" the 
singularity and treating it analytically. 

Constraint-violating instabilities occur because the constraints are not exactly satisfied in a nu- 
merical simulation. Many theorems about general relativity hold only when the constraints 
are satisfied. It is not known what happens when the constraints are not satisfied, and it is 
quite possible that such a constraint violation would grow exponentially. Constraint violat- 
ing modes may also propagate faster than light; they do not have to follow causality. The 
constraint violation can be easily monitored. It is possible to enforce the constraints at any 
time by e.g. solving elliptic equations, but this is expensive. 
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2.6 My goal in this thesis 



Gauge instabilities result from a bad choice of coordinates. With a gauge evolution condition, i.e. 
a choice of lapse and shift, one does not choose the coordinate system directly, but rather 
only its time evolution. Especially when lapse and shift are themselves the result of a time 
evolution equation with its associated boundary conditions, the relationship between them 
and the quality of the coordinate system becomes complicated. Gauge instabilities can theo- 
retically always be cured by a different choice of coordinates, but this is difficult in practice. 
There is no clear-cut distinction between physical and gauge degrees of freedom, which 
makes it difficult to separate between physical and gauge instabilities. 

Instead of classifying instabilities by the affected quantities, one can also try to classify them by 
their sources: 

Instabilities in the formulation are another name for an ill-posed problem. The prototype of an 
ill-posed problem is the backward heat equation u = —u". While the (forward) heat equation 
has a smoothing property, the backward heat equation amplifies any perturbation exponen- 
tially. The amplification time scale depends on the length scale of the perturbation, where 
smaller perturbations are amplified more quickly. 

Instabilities due to boundaries are caused by inappropriate boundary conditions. The advection 
equation li = u' on the unit interval [0; 1] needs a boundary condition at x = 1, but cannot 
have one arbitrarily imposed at x = 0. This follows from the characteristic of the equation. 
However, characteristics for more complicated, nonlinear equations that are not given in a 
first-order form are difficult to find. The characteristics also depend on the formulation and 
the gauge conditions, so that expressing the same physical system with different variables 
or with a different gauge will lead to different characteristics. 

Instabilities in the discretisation are introduced by an unsuitable discretisation scheme. A proto- 
typical example of such an unsuitable scheme for e.g. the advection equation ii = u' is com- 
bining centred differencing with an explicit first order time integrator. Suitable schemes for 
this equation are e.g. using upwind differencing, adding artificial diffusion, using an implicit 
time integrator, or using interpolation instead of finite differencing and time integration. 

2.6 My goal in this thesis 

In this thesis, I am not attempting to find a well-posed formulation for the Einstein equations. I am 
instead using a simple, straightforward formulation that might well not be well-posed (although 
I do not know this). I am adding artificial diffusion to counteract problems caused by this. This 
artificial diffusion cures the discretisation instabilities. 

Also, in this thesis I am not trying to find an especially suited discretisation scheme for the 
Einstein equations. I am using centred differencing and a second order explicit time integrator, 
which are among the most simple tools that can be used for this. First order artificial diffusion 
makes my overall discretisation scheme first order accurate only. 

Instead, in this thesis, I use the first classification of instabilities, and discuss gauge modes. I 
introduce generic gauge fixing conditions which apply in ways equivalent to the constraints. That 
allows gauge modes and gauge instabilities to be identified, and also to be removed by enforcing 
the gauge. 
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3 Gauge conditions 



The ADM evolution equations contain the free quantities lapse a and shift /3', which need to be 
specified prior to a time evolution. Lapse and shift determine how the coordinate system in the 
spacetime is constructed, starting from the current time slice. Such a choice of lapse and shift is 
also commonly called a gauge condition. Unfortimately, the simplest gauge choices will in general 
lead to coordinate singularities, so that more involved conditions are necessary. 

The gauge choice also influences the properties of the system of evolution equations. Hyper- 
bolicity, and the stability properties in general, do crucially depend on the gauge condition. The 
search for good gauge choices forms an important part of the current efforts to obtain stable evo- 
lutions. 

Below, I list several gauge choices that are or were commonly used, and mention some of the 
classifications for gauge conditions that have proven useful in numerical relativity. I discuss the 
difference between gauge evolution conditions and gauge fixing conditions, and point out the proper- 
ties that a good gauge fixing condition needs to have. 

I then describe such a gauge fixing condition that has these properties, and that I will use in the 
following chapters. I also examine the behaviour of its induced coordinate systems and compare 
it to existing gauge conditions. 

3.1 Common gauge conditions 

Clearly the easiest way to specify lapse and shift is by prescribing them as functions of space 
and time, i.e. as functions of the coordinates x' and t. Doing so might be called using an exact 
gauge condition^. Again the simplest possible way to do this is to prescribe a = 1 and = 
everywhere, which is called geodesic slicing and normal coordinates, respectively. For other exact 
gauge conditions, one often uses a known solution of Einstein's equations to obtain expressions 
for the lapse and shift. As mentioned above, such a simple gauge choice usually leads to unstable 
evolutions. The intuitive explanation for these instabilities is that numerical errors accumulate, 
creating gauge modes or a gauge drift, while the gauge condition does not compensate. 

An only slightly more complicated way to specify lapse and shift is by taking the values of the 
primary variables into account. This is often called an analytic gauge condition. In an ADM-like 
system, one specifies lapse and shift as functions of the coordinates x' and t, the three-metric y,^, 
and the extrinsic curvature Kjj. Thus analytic gauges are a superset of exact gauges. The additional 
freedom in choosing lapse and shift can be used to coimteract instabilities. It is also possible to 
solve time evolution equations for lapse and shift. For example, harmonic slicing |BM88| can be 
implemented by setting dtcc = —a^K. Harmonic slicing has singularity avoiding properties similar 
to maximal slicing (see below). 

The most expensive gauge conditions are elliptic gauge conditions, consisting of elliptic equations 
for lapse or shift. They are a subset of analytic gauge conditions. The prototype for this is maximal 

think that this expression was coined by Manuel Tiglio. 
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slicing ISY78I . This gauge condition prescribes K = 0, which leads to the condition Aa = aR for 
the lapse. Other well-known elliptic gauge conditions are minimal strain and minimal distortion 
ISY78I , and the Gamma-freezing shift derived from 3f P' =0 IPolOlI in the BSSN system. 

The gauge conditions currently en vogue in the community are analytic gauges with time evo- 
lution equations. Elliptic conditions have desirable properties, but are considered to be too expen- 
sive. Harmonic slicing and variations thereof can be seen as driver conditions. Instead of solving 
an elliptic equation for lapse or shift, one transforms the elliptic equation into a parabolic or hy- 
perbolic time evolution equation. These evolution equations drive the lapse or shift towards the 
desired value. A parabolic equation acts as a diffusive mechanism that disperses the difference, 
and a hyperbolic equation acts as a (damped) radiative system that transports the difference away. 

Hyperbolic driver gauge c onditions have been successful to some extent for the BSSN system, 
and also for the KST system lKST"'"00l with a densitised lapse. 

As mentioned in the introduction, I do not share the belief that elliptic gauge equations are too 
expensive. While it is obviously true that they are inconveniently expensive, it is often possible to 
optimise them for speed, such as e.g. solving only every n time steps, and using some approxima- 
tion in between. Furthermore, stability is more important than speed. 

3.2 Gauge evolution and gauge fixing 

A gauge condition is in four-space a condition that is applied to the four-metric gyiy. A choice of 
lapse and shift corresponds to a specific choice of goo and gQi and is thus a true gauge condition 
there. 

However, in an ADM-like system, the primary variables are not the four-metric, but rather the 
three-metric and the extrinsic curvature. A gauge condition thus has to be applied to these quan- 
tities. A choice of lapse and shift does not restrict the primary variables, and should therefore not 
quite be called a gauge condition. Because I do not want to break with tradition too much, I call 
it a gauge evolution condition in this text. While a choice of lapse and shift does not eliminate the 
gauge freedom present in the initial data, it does determine how the gauge condition is evolved 
in time. A gauge evolution condition can be obtained from an evolution equation of a primary 
variable, e.g. from dtK = ■ ■ ■ or dtV' = ■ ■ ■. 

A "true" gauge condition is what I call a gauge fixing condition here. A gauge fixing condition is 
a direct restriction on the possible values of the primary variables. This will implicitly also lead to 
a choice for lapse and shift. Maximal slicing (X = 0) is an example of a gauge fixing condition; it 
not only selects a lapse, but also restricts the possible sets of values for the primary variables. 

(Gauge fixing is not a new idea. In fact, gauge fixing is the way in which gauge conditions are 
applied in other fields. Calling a choice of lapse and shift a gauge condition seems to be a relatively 
new convention in general relativity, which was used by Smarr and York in |i SY7$ l|.) 

Most of the gauge conditions mentioned above and used today are only gauge evolution condi- 
tions. This fact does not make them inferior per se, but using them can easily lead to (coordinate) 
singularities, which show up as instabilities. Let me consider, as an example, geodesic slicing 
(a = 1) in a flat spacetime, and let me disregard the shift for now. In Minkowski coordinates, 
geodesic slicing works fine. With different initial data for a flat spacetime, geodesic slicing can 
lead to coordinate singularities, while of course the spacetime is still flat. This is numerically a 
problem when one starts out close to Minkowski coordinates, with the difference being only due 
to numerical errors. 
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Unexpected gauge singularities, i.e. gauge instabilities, cannot happen with a gauge fixing con- 
dition. A gauge fixing condition imposes a condition on a primary variable, e.g. specifies the 
values for K, and then selects a lapse so that this condition continues to hold. Because K is speci- 
fied, it is impossible to accidentally select a lapse that makes K diverge, as in the example above. 
The only way to produce a gauge singularity is by specifying a diverging K to begin with. 

This does not prevent accumulated numerical errors, which could still make K diverge. But the 
gauge fixing condition can be monitored during the evolution, so that one can judge the quality 
of the simulation. This is equivalent to the way one handles the constraints. Furthermore, a gauge 
fixing condition can also be directly enforced when numerical errors have accumulated, eliminat- 
ing gauge instabilities altogether. 

The interplay between lapse and shift and gauge singularities is further complicated by the 
boundary conditions. One needs boundary conditions for the primary variables, and also for 
lapse and shift if they are evolved in time. Without a proper gauge fixing condition to start from 
it is prohibitively difficult to ensure that a certain choice of lapse and shift does not, under the 
influence of numerical errors, lead to a gauge singularity at the boundary. With a gauge fixing 
condition, such instabilities are not possible, or can at least be easily detected. 

The only gauge fixing condition that is in use today in nonlinear three-dimensional numerical 
relativity is maximal slicing. 

A close candidate is AEI's Gamma freezing, which sets 3f F' = and derives an elliptic shift con- 
dition from that. Placing a condition on a time derivative is not directly a gauge fixing condition, 
which requires placing a condition on the variable itself. But setting a time derivative to zero can 
be viewed as fixing the variable at its initial value, and I therefore want to count it as gauge fixing. 
However, Gamma freezing does not try to keep the value of F' consistent with y,y, i.e. the new 
constraints arising from introducing F' as additional primary variable IBS99I are not controlled. 
(These constraints seem to be important, because they seem to be growing without bound when 
they are not controlled.) Re-calculating F' from y,^- destroys the gauge fixing properties of Gamma 
freezing. 

Minimal distortion is not a gauge fixing condition. This condition minimises the distortion, 
which is the change of shape during a time step. It does not try to attain any specific coordinate 
shape (e.g. for a sphere), and is therefore not a gauge fixing condition. It restricts only the time 
derivative of the three-metric, and not the three-metric itself. 

One interesting case of a gauge fixing condition is area locking. Area locking was introduced in 
spherically symmetric situations, and it leads to very stable evolutions in this case. It sets dtTee = 
0, which leads to a condition for the shift component 13''. As above, I count this as gauge fixing. 
Unfortunately, area locking was never successfully generalised to three dimensions. 

3.3 Good gauge fixing conditions 

A gauge fixing condition is a condition that acts directly on the primary variables, i.e. the three- 
metric and the extrinsic curvature for an ADM-like system. It has thus properties similar to the 
constraints. Gauge fixing conditions are a superset of gauge evolution conditions which select a 
lapse and a shift only. 

A good gauge fixing condition has to have several properties: 

1 . Of course, it has to leave enough degrees of freedom so as to not constrain any physics, even 
in the full nonlinear case. No physical degree of freedom must be removed from the system. 
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2. It has to lead to a unique prescription for lapse and shift, so that the gauge condition will 
also hold for the neighbouring time slices. 

3. It should be easy to evaluate, so that one knows whether the gauge condition still holds, or 
by how far one has left this gauge. 

4. There should be a reasonable way to enforce this gauge, so that one can counteract numerical 
errors that have accumulated over time, or create initial data in this gauge. 

5. It should be independent of the constraints. That is, the condition should be independent of 
the variables that are changed when the constaints are enforced, so that one can at the same 
time enforce the gauge condition and the constraints. 

6. It should have an understandable physical interpretation, so that the physicist can apply his 
or her intuition. 

7. It should be general enough so that many kinds of initial data can be prepared with this 
gauge, and many spacetimes can be described in it. 

I do count solving elliptic equations as being "reasonable" in the context of the above. There 
are kinds of equations that are much more unwieldy than elliptic ones. And there exist efficient 
methods for solving elliptic equations with a cost that grows only linearly with the problem size. 
I do not strive for the fastest possible implementation here, I only want something "reasonable". 

Maximal slicing almost satisfies these properties, but it is not a complete gauge choice. It selects 
a slicing, but does not determine the three-coordinate system. Additionally, the condition X = is 
too constraining in practice — it does e.g. not permit a horizon-penetrating static slicing of a static 
black hole (see e.g. LCoo02. section 3.1.1]). 

3.4 My gauge fixing condition 

My goal in this thesis is to present a gauge fixing condition that can be used to simulate black 
hole spacetimes. This condition has to have the properties that I list above. As mentioned in 
the introduction, I also want a simple system. The currently promising systems (e.g. the BSSN 
variants) all contain many subtle rules that have to be followed closely. The gauge condition 
should be conceptually as simple as possible, even if it is expensive to use. 

For my gauge condition, I choose as gauge variables the trace K of the extrinsic curvature and 
certain combinations of partial derivatives of the conformal three-metric, which I call F; after Shi- 
bata and Nakamura iSN95l- My gauge condition then consists of prescribing, in advance, values 
for these quantities as functions of space and time, i.e. depending on x' and t only. 

Prescribing the trace of the extrinsic curvature leads to an elliptic equation for the lapse (see 
appendix lA.2.2> . Given that, it is clear that a choice of K has all of the above properties. Interpreting 
i<C as a gauge variable is a common choice. It is a straightforward generalisation of maximal slicing 
ISY78I ■ and I think that this gauge variable needs no further justification. 

My gauge condition acting on the metric is a bit more involved. I first define (for technical 
reasons) the traceless conformal three-metric hjj as 

% = 7ij-l^rj^"fkl ■ (3.1) 
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The bar ■ indicates here that the indices are to be raised and lowered with the coordinate metric 
5jj. My gauge variable F, is then defined by 

Fi = h,j.j ■ (3.2) 

That is, my gauge quantity F,- is closely related to the divergence of the conformal three-metric yij- 
The fact that 1 use the traceless part only has only technical reasons, having to do with the method 
used to enforce the gauge condition, which is explained further down. 

This gauge condition is similar to the transverse-traceless gauge used in linearised general rel- 
ativity, where one often requires of the metric perturbation the condition hjj^j = 0. (Here /z/y is 
a three-metric perturbation, which is different from the variable S/y above.) However, 1 do not 
require that F, = 0, but only that it be prescribed in advance. 

This gauge variable is also very similar to, but slightly different from the auxiliary variable F, 
introduced by Shibata and Nakamura ISN95I as Fj = y^y y. The BSSN system has instead P' = 
y^^ Fyj. IBS99I . which is also closely related to the Shibata-Nakamura F, . 

It is by design that the gauge condition is not covariant. The gauge defines the coordinate 
system, and the use of a specific coordinate system is what makes a calculation be not covariant. 
That means that it does not make sense to speak of a covariant gauge condition — such a gauge 
condition could not single out a specific coordinate system. 

The choice of the gauge variable F, constrains the possible shapes of the three-coordinate sys- 
tem. The gauge condition judges the quality of the current coordinate system, and prevents a 
deterioration by placing the grid points appropriately. 

As described in section \TM the metric components y„ (no summation implied) describe the 
physical distance between grid points in the i direction, while the components y,j (with i j) 
describe the physical angle between the i and ; coordinate directions. The direct metric divergence 
y/j y hence describes the spatial change of the behaviour of the grid lines in the i direction — the 
change in length, plus the change in the angles with the other grid lines. 

The physical meaning of the components of the conformal traceless metric hjj is more complex, 
but will be similar in nature. Hence enforcing a certain hij j means to enforce a certain spatial rate of 
change of the coordinate system. Enforcing the gauge changes the metric, and therefore effectively 
moves the grid points until they have the "right" location with respect to their neighbours. 

As one can set the quantity F; to an arbitrary function, this gauge choice can be used for all 
initial data. Fixing F, leads, via the time evolution equation of y,j, to an elliptic equation defining 
the shift (see appendix lA. 2. 7L so that using this shift preserves the gauge condition. That means 
that this gauge condition is generic, can be used for all spacetimes, and does not constrain any 
physical degrees of freedom.^ 

The resulting elliptic equations for the lapse oc and the shift 13' also require boundary conditions. 
These boundary conditions are part of the gauge condition. One can use the boundary conditions 
for the lapse to slow down or speed up the time evolution, or to tilt the time slice. The boundary 
conditions for the shift can be used to move the boundary further in or out, or to move or rotate the 
whole simulation domain. This can be used to set up co-rotating coordinate systems, for example. 

By introducing a vector potential W, for the traceless conformal metric hjj, one can enforce this 
gauge condition in a manner equivalent to enforcing the momentum constraint on the extrinsic 

^In this I assume that this elliptic equation has a solution. This depends e.g. also on the boundary conditions. For the 
cases 1 am interested in, one can assume that the domain is simply connected and has Dirichlet boundaries. I have not 
proven that a solution exists, but numerical tests support my assumption. 
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curvature.'' One can reconstruct the conformal metric y,y from its traceless part through the con- 
dition dety,y = 1, which determines its trace. Furthermore, as enforcing the constraints in the 
York-Lichnerowicz formahsm (see section l4!2l does not change the conformal metric or the trace 
of the extrinsic curvature, it preserves this gauge conditions. 

This means that this gauge has all the properties (1) - (7) above for a good gauge fixing condition. 
I will use it in the following chapters and examine its numerical properties. 

Another, potentially more elegant gauge condition on the metric could be F, = fijj- Another, po- 
tentially more elegant way of enforcing the gauge condition could be to apply a four-dimensional 
coordinate transformation. A coordinate transformation would have the advantage that the con- 
straints are not influenced. I did not investigate these alternatives. 



^The existence and uniqueness properties of ttie York-Lichnerowicz method thus are trivially valid for this method as 
well. 



26 



4 The system of equations 



4.1 The variables 

I chose not to use the standard ADM system to represent a spaceUke hypersurface, but rather a 
variant of a conformal traceless ADM (ctADM) system. Such a system can more easily be used to 
enforce the constraints, and furthermore it also already makes the gauge degrees of freedom more 
explicit. 

My primary variables are 

• the conformal factor xp, with i/)^^ = det y 

• the conformal three-metric fjj, with t/^^7,y = 

• the trace K of the extrinsic curvature, with K = y'^Kij 

• the traceless conformal extrinsic curvature with ip^^Aij = Aij. Here A/y is the traceless 
extrinsic curvature, i.e. Ajj = Kjj — ^YijK 

• the divergence F, of the traceless conformal three-metric, with F, = hjj j. Here hjj is the 
traceless conformal three-metric, i.e. hjj = fjj — ^SijS^^f^-i 

In the above, y, j and j are the three-metric and extrinsic curvature of the standard ADM system. 
This system implicitly fulfils the constraints dety,^- = 1 and trace A/y = y'^A/^ = 0. I use a tilde ~ 
for quantities that are connected to the conformal metric y,y rather than the physical metric y,y, and 
a bar " for quantities that are connected to the coordinate metric The time evolution equations 
for this system are derived in appendix lA.2.21 the constraints in these variables in appendix IA.2. 31 
This system, which I call the TGR system, is similar to the BSSN INOK87I In089 '. 'SNgS . "55991 
system. The main differences are that it has the conformal factor i/; instead of its logarithm (p = 
log xp as primary variable, and that the traceless cortformal extrinsic curvature is scaled with a 
different power of the conformal factor; while the TGR system uses the scaling factor the 
BSSN system uses the same factor ip'^ as for the conformal metric. It also has F, instead of F', which 
is similar but different. 

Compared to the original ADM system, the TGR system has the advantage that it can more 
easily be used to enforce the constraints. Enforcing the constraints using the York-Lichnerowicz 
formalism IYor73l IVor74l ICoo02l section 2.2.1] requires a decomposition into a system with a 
conformal metric and a traceless conformal extrinsic curvature. With the standard ADM system, 
this requires variable transformations before and after each enforcing, which leads to additional 
discretisation errors. In fact, the very system used by Cook in ICoo02l section 2.2.1] motivated 
my choice of variables. The BSSN system would also suitable for solving the constraint equations, 
although only with a slightly different method. 
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4.2 Enforcing the gauge and the constraints 

In order to enforce the constraints with the York-Lichnerowicz method, one introduces a vector 
potential V, for the traceless conformal extrinsic curvature. This transforms the constraint equa- 
tions into a set of coupled nonlinear elliptic equations for the conformal factor and this vector 
potential, which can then be solved. This is described in appendix IA.2.41 The gauge condition 
can be enforced in a similar way by introducing a vector potential W, for the traceless conformal 
metric. This is described in appendix lA.2.61 

With given initial data, boundary, and gauge conditions, the time evolution of a physical system 
is defined by the time evolution equations for the above system. The constraint and the gauge 
equations can also be used to enforce the contraints and the gauge, i.e. to create initial data, or to 
counteract numerical errors during the time evolution. 

Enforcing the constraints and the gauge conditions has to happen in a certain order. This or- 
der is necessary, because e.g. enforcing the gauge condition for F, on y,y changes the constraints. 
Therefore, the gauge condition for F, has to be enforced before the constraints are enforced. Lapse 
and shift are calculated after all enforcing. One possible order is the following: 

1. Enforce gauge, part I: set K = K* 

2. dety,-^-: Rescale fjj such that dety,j = 1 

3. Enforce gauge, part II: change jjj such that F, = F*, preserving det y,y 

4. trace A,y: Change such that trace A/^ = (depends on fjj) 

5. Enforce constraints: change xp and A,^ such that H = and M, = (depends on fjj, K, and 
Aij) 

6. Calculate lapse a and shift (depends on xp, fjj, K, and A,y) 

In the above, values marked with an asterisk ■* are values prescribed by the gauge condition. 

In principle, enforcing the gauge and the constraints is not necessary for a perfect time evolution. 
In practice, it is necessary or advantageous to counteract numerical errors. It might be possible to 
change these elliptic equations into hyperbolic ones which would be much cheaper to solve; I have 
not tried this. 
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Boundary conditions are necessary while integrating the evolution equations in time, and for solv- 
ing the elliptic gauge and constraint equations. They are needed at the outer boundary, and if 
excision (see below) is used, also at the excision boundary. 

5.1 Outer boundary 

Most astrophysically interesting spacetimes are asymptotically flat IWal84l chapter 11]. That 
means that at large distances from the origin, the spacetime is flat plus some perturbation that 
falls off at least with 1 /r. Of the spacetimes considered in this thesis, only some numerical test 
problems are not asymptotically flat. 

For single black hole and binary black hole collision simulations, the outer boundary should 
ideally be located at spatial or null infinity. In that case, no gravitational radiation will enter or 
leave the simulation domain, which makes the boundary condition easy to implement. 

5.1.1 Location of the outer boundary 

In a numerical simulation, the outer boundary will be at a finite location (and not at spatial infin- 
ity), where the four-metric is not flat.^ It can usually be assumed to consist of a known backgroimd 
metric that is superposed with small perturbations. These perturbations can be gravitational 
waves, or can be gauge modes, or can be constraint violating modes that arise from numerical 
errors. 

It is theoretically possible to put the outer boundary indeed at spatial infinity by choosing a 
suitable coordinate system. However, one would have to get there with a finite number of grid 
points, and the gravitational waves would be backscattered by the change in resolution. One 
would have to make sure that the gravitational waves have been absorbed or dissipated before 
they reach a region where they cannot be resolved any more. Given this, there seems to be no real 
advantage in putting the outer boundary at spatial infinity. 

5.1.2 Physical boundary conditions 

Astrophysically, one wants an outgoing radiation boundary condition at the outer boundary. All 
gravitational radiation that reaches the outer boundary should be let out, and no gravitational 
radiation should be sent in or reflected back. This means that the outer boundary has to be in a 
region where one can meaningfully speak of gravitational radiation "on top of something else", 
i.e. one needs to be close to a known background solution, and the radiation ampHtude has to be 
small. 

^It is possible to put the outer boundary at null infinity or beyond by using a conformal approach IFra98l . I will not 
consider this here. 
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A true radiative boundary condition would require splitting the quantities at the boundary into 
a background solution containing the total mass and spin, plus a linearised perturbation repre- 
senting the gravitational waves, separated into their incoming and outgoing modes. This would 
be a rather complex and time consuming process. People have instead contemplated matching a 
perturbative formulation to the nonlinear formulation at the outer boundary | Win02|. This should 
in principle lead to a very clean boimdary condition, but is hampered in practice by the coordinate 
transformations necessary to match the formulations. 1 know of no stable implementation using 
this approach. 

One approximation is to use some kind of sponge layer near the outer boundary. Within this 
layer, one introduces an increasing absorption coefficient. It emits no radiation by itself, absorbs 
the waves reaching the layer, and reflects almost nothing back. On the outer side of t his sponge 
layer one imposes a Dirichlet boundary condition. Perm State uses a blending zone IBCG+OOI . 
which is similar to this approach. 

One other approximation is to assume that the waves are spherical about the origin. AEl uses 
conditions like this. Such radiative boimdary conditions assume that the evolved quantity is an 
outgoing radial wave, satisfying f{t, r) = foo + u{r — t)/r, where the constant /oo is given, and the 
fimction u is determined implicitly at each boundary point. 

For elliptic equations, the usual boimdary condition is a Robin boundary condition, which en- 
forces a certain falloff towards infinity, i.e. a condition f{r) = foo + C/r" with given /oo and n. 
Here foo is the desired value at infinity, and n is the falloff power. The constant C is determined 
implicitly at each boundary point. 

If the solution at the outer boundary is known, which is usually not the case, then one can use 
Dirichlet boundary conditions. For waves, a constant-in-time Dirichlet boundary is a reflective 
boundary, so that waves cannot leave the simulation domain. This makes Dirichlet boundaries 
unsuitable for realistic applications, but also well suited for test problems, e.g. to study the stability 
of a formulation. 

5.1.3 Periodicity 

It is, especially for test problems, sometimes convenient to assume periodicity. Periodicity means 
that there are no real boundary conditions; this allows one to examine the system of evolution 
equations independent of any problems that might be caused by the boundary conditions. For 
example, the first two stages of Winicour's robust stability test |SGB WOO SSW02 1 involve periodic 
boundaries. Additionally, periodicity allows one to more easily examine the long-term behaviour 
of dynamic systems, because the interesting features carmot leave the simulation domain. 

Unfortunately, periodicity poses a severe problem for elliptic equations. Their solution is deter- 
mined by the boundary condition, and without a boundary their solution is not unique any more, 
or might cease to exist. If the solution is not unique any more, then one can impose a pseudo 
"boundary condition" at a single point in the domain to select a solution. 

The elliptic equations arising in the TGR system seem in fact to be ill-posed when given periodic 
boundaries. That is, they do not admit a solution any more. This can already be seen at the 
example of the Poisson equation Acp = p with periodic boundaries. For p = 0, one needs to 
define <p at one point to make the solution unique. For p 7^ 0, the shape of is a parabola, 
and is therefore clearly not periodic. Similarly, the Hamiltonian constraint equation in the York- 
Lichnerowicz formalism with periodic boundaries seems to be ill-posed when Kjj is not zero. For 
Kjj = 0, a solution exists, but is not unique. 
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This prevents many interesting test cases from being run while enforcing the gauge conditions 
and constraints. It is in some cases possible to convert these test cases into almost equivalent ones 
with Dirichlet boundaries. A travelling wave with periodic boundaries can e.g. be changed into a 
standing wave with periodic boundaries, and that system is then very similar to a standing wave 
with Dirichlet boundaries at the wave knots. And fortunately, test cases with Dirichlet boundaries 
are often considered to be more difficult to pass than those with periodic boundaries, i.e. they form 
stronger test cases. 



5.1.4 Gauge and constraint boundary conditions 

A physical boundary condition does not define how to handle the gauge and constraint violating 
degrees of freedom. Physically, the constraints have to be satisfied. Numerically, there will always 
be errors. One has to try to reduce the violation of the constraints in some way, and this mechanism 
will interact with the boundary condition. 

Similarly, after fixing the gauge, the gauge will numerically only be satisfied approximately. The 
boundary conditions have to deal with violations of the gauge and the constraints, meaning that 
there also has to be a boundary condition for the gauge and constraint violations. Unfortunately it 
is not clearly defined just which (combination) of the evolved or solved-f or quantites contain these 
gauge and constraint violating degrees of freedom. 

In my system, the variables ip, fjj, and Ajj, are evolved in time, while K and F, are prescribed in 
advance by the gauge condition. It is not necessary to evolve the conformal factor ip in time, as it 
can also be calculated from the Hamiltonian constraint. I do evolve it only to have a better initial 
guess for enforcing the Hamiltonian constraint, and in order to have more freedom in choosing 
the boundary conditions. Depending on the problem, I use either Dirichlet or radiative boundary 
conditions for the time evolution of these evolved variables. 

It is clear that a Dirichlet boundary condition, even if it is consistent with the equations, is not 
suitable for an astrophysical problem. The radiation that is reflected at the outer boundary makes 
the result unphysical. But as such reflected radiation is considered to introduce instablities, a 
Dirichlet boundary is actually a stronger stability test. I would use a sponge layer or a radiative 
boundary condition if I were to run a simulation with a domain large enough and a resolution fine 
enough to permit meaningful extraction of gravitational radiation. 

The variables W,, Tp, and Vj need boundary conditions when the gauge and the constraints 
are enforced by solving elliptic equations. I use either Dirichlet or Robin boundary conditions 
for them. I impose W, = and V, = on the vector potentials. For the conformal factor xp I 
distinguish between initial data and time evolution. For initial data I use either a Dirichlet or 
a Robin boundary condition. During time evolution, I use a Dirichlet or a radiative boundary 
condition. 

I also apply Dirichlet or Robin boundary conditions to the lapse a and the shift p' . When using 
a Dirichlet boundary condition, I retain the boundary values from the initial data. 

The lapse and the shift determine the coordinate system that is constructed for the spacetime 
during time evolution, and therefore their boundary values are important and have a clear mean- 
ing. Their boundary condition is part of the gauge condition, as described in section lT^ 
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5.2 Excision boundary 

Black hole spacetimes contain singularities. These singularities can, of course, not be directly 
treated numerically. Several methods have been used to treat them in a 3 + 1 split. 

The likely oldest and mathematically most elegant ansatz is maximal slicing ISY78I . For this, one 
chooses the coordinate time in such a way that it reaches the singularity only for f — > cxd. The 
relation between the coordinate time t and the proper time t of an observer can be chosen almost 
arbitrarily through a suitable lapse function. With maximal slicing it is nevertheless the case that 
the whole spacetime is eventually covered by the time coordinate. 

This elegant ansatz has the disadvantage that the time coordinate in the neighbourhood of the 
singularity proceeds more and more slowly, which eventually leads to a distortion of the coordi- 
nate system ("grid stretching"). Although the coordinate system stays well defined all the time, 
certain metric components start to grow without bound. This can in the end not be treated nu- 
merically any more. Maximal slicing can thus describe a black hole only for a certain amount of 
coordinate time before a code crashes for numerical reasons. There are also different slicings, such 
as harmonic |BM88J or 1 -|- log LACM''~95J slicings, with similar singularity-avoiding properties, 
and which also have similar disadvantages. 

A different ansatz was developed at the AEl and is called punctures fBB97l. At the singularity, 
certain metric components are infinite. This behaviour can be described by decomposing the three- 
metric Yij into a conformal factor i^i and a conformal three-metric 7,^. For certain initial data (such 
as e.g. Brill-Lindquist black holes; see section I6.3.1> . the conformal three-metric remains finite 
over the singularity, and only the conformal factor diverges. The conformal factor for the initial 
data is known analytically, and by certain gauge choices (maximal slicing (X = 0) and normal 
coordinates (/3' = 0) at the punctures), the time derivative of the conformal factor can be set to 
zero. Additionally, one chooses a numerical grid such that the singularity does not coincide with a 
grid point. Through that, the code never has to calculate any diverging quantities. The conformal 
factor and its spatial derivatives are known analytically and can be provided with arbitrarily high 
accuracy. 

However, the restriction on the gauge conditions has certain disadvantages; e.g. the position 
of the singularity cannot change in time. When two black holes coalesce, their physical distance 
decreases, while the coordinate distance stays constant with this method. This leads to similar 
problems as with maximal slicing. But a single black hole, which is an important test problem in 
numerical relativity, can in principle be described without problems with pimctures. 

Another way to treat singularities is excision. It consists of removing a certain region around the 
singularity from the simulation domain. Because there is, under certain reasonable assumptions, 
always an event horizon enclosing a singularity, this excising is possible without influencing the 
part of the domain outside the horizon. In the language of hydrodynamics, the excision or inner 
boundary is a supersonic outflow boundary for the physical degrees of freedom. Outflow means 
here "out of the simulation domain", i.e. in the direction towards the singularity, further into the 
black hole. 

In spherical coordinates, excision is quite customary, and it poses no problems, as the boundary 
can be put at an r = const coordinate plane. Unfortunately, spherical coordinates are not suitable 
for describing two coalescing black holes, which is why one usually uses Cartesian coordinates for 
that. The excised region on a Cartesian grid can not be spherical any more, but will rather look 
like a sphere that has been built from Lego blocks (see figure l5Ai . This irregular shape leads to 
problems when implementing the excision boundary condition. 
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Figure 5.1: Left: Schematic picture of an excision region. The circle depicts the apparent hori- 
zon. The grey grid cells are excised. The excision boundary has to follow the grid cell 
boxmdaries. Right: A cut through the equatorial plane for an M = 1, fl = 1 black 
hole in Kerr-Schild coordinates. There is not much coordinate space between the event 
horizon and the ring-shaped singularity. 



It is generally assumed that excision is the most promising way to treat singularities, and that it 
will be widely used as soon as good boimdary conditions are known for the excision region. 

I will only consider excision boundaries in the remainder of this chapter. I will not consider 
pimctures at all, and in order to use maximal slicing, one needs of course no special treatment of 
the singularity. 

5.2.1 Location and shape of the excision boundary 

In order to excise a region of space, one generally has to find a closed two-surface that is guar- 
anteed to be both within the event horizon, and away from the singularity. One surface that is 
guaranteed to be within or on an event horizon is an apparent horizon (see section I7!6l . For nu- 
merical reasons, one wants to stay several layers of grid points away from the event horizon, and 
also sufficiently far (as far as possible) from the singularity. It is thus necessary to locate apparent 
horizons to define or at least check the location of the excision boundary. 

Choosing the right shape for the boundary is another rather difficult issue. Spherical or elliptic 
shapes are difficult to handle, because the do not have a normal direction that is aligned with the 
grid. Such normal directions are necessary for extrapolations. On the other hand, using a cubic 
boundary is not always possible. For example, the singularity of a rotating black hole in Kerr- 
Schild coordinates (see section l6. 2. 1> is ring-shaped, and the space available between the ring- 
shaped singularity and the elliptic horizon shape is too small to permit using a cubic boundary (see 
figure lSTl . One thus has to deal with boundary shapes that are not aligned with the computational 
grid. Some possible shapes include boxes "with the corners cut off", which are polyhedra with 18 
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or 26 faces.^ 

5.2.2 Boundary conditions 

Excision boundaries are very different from outer boundaries. Apart from the fact that they are 
more difficult to handle numerically because of their irregular shape, they are located in the strong 
field regime, and hence one cannot assume that any approximation holds. Luckily, there is an 
event horizon present, and the fact that no physical information can escape from the horizon 
means that one has some freedom in managing this boundary. As long as one stays within the 
limits of physics, one can choose any boundary condition without influencing the physics of the 
spacetime outside the horizon. 

Yet there are things that do escape from the horizon. Information about the gauge and about 
constraint violations does not carry physical information. It can propagate faster than light, and 
can also leave the horizon. There is a priori no reason to assume that it would not do so. Care has 
to be taken in separating the individual quantities, which is much more difficult than at the outer 
boundary, where everything can be assumed to be a small perturbation of a known background. 
The inner boundary is a place where a hyperbolic formulation would definitively be of help. 

While solving elliptic equations, it is consistent to use (arbitrary, constant in time) Dirichlet 
boundary conditions for the quantities V„ W„ a, and j3\ Here the boundary values of the shift 
have to prevent the grid boundary from falling further into the black hole and getting too close 
to the singularity, or from rising out of it and leaving the event horizon (see also section l23l . The 
latter would invalidate the assumption underlying excision and must therefore not be permitted. 

While integrating in time, the boundary value for the conformal factor ip has to be determined 
by a Dirichlet boundary condition from a known solution, or by extrapolation. The quantities y,y 
and Ajj carry a mixture of physical, gauge, and constraint information. The physical information 
has to be advected out of the simulation domain. Using a Dirichlet boundary condition for the 
physical degrees of freedom is consistent only if one uses the correct values. Of course, the correct 
values are only known if one tests the code with an analytic solution, so this is not possible for 
real- world cases. 

A supersonic outflow boundary condition can be realised by extrapolation. Unfortunately, ex- 
trapolating the conformal metric causes problems, because the extrapolated conformal metric will 
normally not satisfy the constraint det y/y = 1 any more. It is not a good idea to enforce this con- 
straint by rescaling the whole metric, because this will change certain gauge degrees of freedom in 
unpredictable ways. Allowing for arbitrary gauge changes at the excision boundary might allow 
gauge modes to appear. 

Instead, the gauge degrees of freedom have to be fixed, according to the plan I set out in sec- 
tion |^| The obvious gauge degree of freedom to fix at the inner boundary is the location of the 
boundary itself. That is, I want to clearly define how the inner boundary moves inward or out- 
ward (again, see also section Therefore I decompose the metric components fjj into those 
parallel to and those normal to the boundary. I enforce the constraint det y, j = 1 by rescaling only 
those components that are normal to the boundary. This is also consistent with the fact that gravi- 
tational waves cause transverse length changes only, because these lengths will not be changed by 
the rescaling. 

The traceless conformal extrinsic curvature happens to be less critical. It does not contain any 
gauge degrees of freedom. By extrapolating Ajj one violates the constraint A. = 0, especially 

think this was first suggested by Miguel Alcubierre in 2000. 
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because the conformal metric used to calculate that trace is also extrapolated, and hence changes 
as well. Both the extrapolation of A,y and enforcing A' = later on will change the boundary 
condition for the constraints. However, the constraints can be enforced no matter what the bound- 
ary values. Numerical experiments show that the handling of the traceless conformal extrinsic 
curvature is not critical and does not lead to instabilities. 

I assume that the constraint violation error that is introduced by extrapolating A, j is less malev- 
olent than the gauge error introduced by extrapolating fij. This gauge error has the ability to 
change the extent of the simulation domain, and that is not acceptable. On the other hand, the 
constraint error could lead to certain components of Aij growing without bound, although the 
constraints would stay satisfied. Luckily, that seems to not happen. 

With the the above boundary conditions for the excision boundary, the time evolution of the 
TGR system is consistent, and all degrees of freedom are fixed. 
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6 Initial data 



It is a nontrivial task to generate initial data for a black hole simulation. There are many analytic 
solutions for single black holes, but the proposed methods to construct spacetimes with two or 
more black holes either restrict the possible configurations, or require solving elliptic equations. 
Multiple black hole initial data usually need to be interpreted in terms of their ADM mass and 
spin fWei75l chapter 6.6], and the apparent horizons (see chapter]^) that can be found. 

Initial data are usually presented in the ADM variables y,^- and even if they are actually cal- 
culated in other variables. This has the advantage that the form of the initial data is independent 
of the formulation used in the evolution, and facilitates exchanging initial data between different 
formulations. The transformation between the ADM and other variables is usually straightfor- 
ward. 

I present three major kinds of initial data in this chapter. I start with data without black holes. 
These data often have an analytic solution for all times and can thus be easily used as test cases. 
I then present various coordinate systems for static or stationary single black holes. They are also 
mostly test cases of little astrophysical relevance. Last I present initial data for multiple black 
holes, for which no closed form solution exists for later times. 

I follow the convention of giving names not to spacetimes, but rather to coordinate systems. This 
is motivated by the fact that different coordinate systems behave very differently numerically. On 
the other hand, the fact that two different coordinate systems might describe the same physical 
spacetime is largely irrelevant for me, because it cannot easily be tested. I also hope that this helps 
to prevent confusion. I thus distinguish between flat space (a spacetime) and a Minkowski metric 
(one possible coordinate system for it). Similarly, there are static black holes, of which Schwarz- 
schild is just one possible metric, with zero-spin Kerr-Schild and Painleve-GuUstand being others. 

6.1 Data without black holes 

6.1.1 Minkowski metric (flat space) 

The Minkowski metric is the simplest case. It describes flat space in Cartesian coordinates. It has 
gfiv = r\^y, or y,y = bip Kij = 0, a = 1, and 13' = 0. 

6.1.2 Weak Bondi wave (linear planar wave) 

A weak Bondi wave is a planar gravitational wave with a small amplitude. This is a solution to 
the linearised Einstein equations, presented e.g. in f MTW73l chapter 35.9, eqn. (35.32)]. Assuming 
that the wave propagates in the z direction, the ADM quantities are given by 



ryy 



Txx 



= 1 + b 
= 1-b 



(6.1) 
(6.2) 
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Tzz = 1 (6.3) 
Kxx = -\h (6.4) 

= \b (6.5) 

Xzz = (6.6) 

where all other components of jij and Kjj are zero, and a = 1 and (3' = 0. The free parameter 
function b{t, z) can e.g. be chosen as 

b = Asinfc(z-t) (6.7) 

h = dtb = -kAcosk{z-t) (6.8) 

for a wave propagating in the z direction, or as 

b = A sinfcz cosfcf (6.9) 
b = dtb = -kAsinkzsinkt (6.10) 

for a standing wave, where in all cases 1^1 <C 1 has to hold. The parameters A and k determine 
the amplitude and wave number of the weak Bondi wave. This metric is a solution of the Einstein 
equations only for \A\ <C 1 . 

6.1.3 Gauge pulse (nonlinear planar gauge wave) 

A gauge pulse is a planar nonlinear gauge wave. This is a solution to the full vacuum Einstein 
equations which does not contain gravitational waves. It represents a flat spacetime in a strange 
coordinate system. Assuming that the wave propagates in the z direction, the ADM quantities are 
given by 

Yx. = 1 (6.11) 

Tyy = 1 (6.12) 

Tzz = exp{b} (6.13) 
1 (b 



Kzz = --exp^-\ b (6.14) 



where all other components of y,y, and /3' are zero. The free parameter function b{t, z) can e.g. 
be chosen as 

b = Asinfc(z-f) (6.16) 
b = dtb = -kAcosk{z-t) (6.17) 

for a wave propagating in the z direction, or as 

b = A sinfcz coskt (6.18) 
b = dtb = -kAsinkzsinkt (6.19) 

for a standing wave. The parameters A and k determine the amplitude and wave number of the 
gauge pulse. The amplitude A can be arbitrarily large. 
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6.1.4 Brill wave 

A Brill wave is a nonlinear, axially sjnnmetric gravitational wave IBri59l . A Brill wave can be 
strong enough to form a black hole. This is an interesting test case, insofar as a black hole can form 
from scratch, without a singularity or matter being present in the beginning. 

In contrast to the solutions presented above, which are valid for all times, the Brill wave metric 
given below describes only initial data at t = 0. The metric can be written as 

ds^ = xp'^ [e^'?(dp2 +dz^)+ p^d(pA (6.20) 

with the free function q{p,z) than can e.g. be chosen as 

q = Ap^e^'^ (6.21) 

with the amplitude A. p is the cylindrical radial coordinate with = + y^. 

These initial data are then chosen to be time-symmetric with 9ty/y = 0. Together with the gauge 
choice /3' = this leads to the extrinsic curvature K;^ = 0, satisfying the momentum constraint 
identically. The conformal factor t/j has to be chosen (solved for) so that the Hamiltonian constraint 
is satisfied. One usually uses a Robin type boundary condition for this. 

One usually uses maximal slicing {K = 0) and normal coordinates (/3' = 0) when simulating a 
Brill wave. One also usually uses a Robin boundary condition for the lapse. 



6.2 Single black hole data 

There are many coordinate systems that describe single black holes. I chose to test my formula- 
tion with three, namely Kerr-Schild coordinates, Painleve-Gullstrand coordinates, and harmonic 
coordinates. The original Schwarzschild coordinates are not suitable here, because they do not 
penetrate the horizon. That would not allow me to use an excision boimdary condition (see sec- 
tion l5.2> , which has to be applied inside the event horizon. 



6.2.1 Kerr-Schild coordinates 

For single black hole test runs I prefer the Kerr-Schild coordinates^ . This is a static (or stationary, 
when there is a nonzero spin) solution with a spin that can be specified freely. The slicing does 
intersect the singularity, and the singularity in the time slice is a point for a non-spinning, and is 
ring-shaped for a spinning black hole. The singularity diverges with only, which still allows 
reasonable numerical resolutions to be used. 

Kerr-Schild coordinates |Coo02 section 3.3.1] are usually expressed in an elliptic coordinate 
system, where the r coordinate is given by 

p2 = r^ + a^[l-pj (6.22) 

^For some historic reason, these coordinates are sometimes also called "ingoing Eddington-Finkelstein coordinates" in 
the community, although the real iEF coordinates are different. 
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with the usual radial coordinate p, i.e. = x'^ + + z^. The four-metric can be written as 



gn-v = ri^^ + 2Hl^ly (6.23) 
where t]^^ is the Minkowski metric, H is a scalar, and is a null vector, which are given by 



H 



rx + ay 

r^ + a^ 
ry — ax 

r^ + a^ 

z 

r 



(6.24) 
(6.25) 
(6.26) 

(6.27) 

(6.28) 



where M is the black hole mass, and a is its spin about the z axis. The three-metric, lapse, and shift 
then follow as 



Yrr 

Trip 

Tee 

a 
13' 



1 + 



2Mr 



2Mr 



J J 2Mr , . , 
r +a A ^a sm i 



sin^e 



2Mr 



,2Mr 



(6.29) 

(6.30) 

(6.31) 
(6.32) 

(6.33) 
(6.34) 



The event horizon is located at r = M -|- V Nf- — a^. 

The extrinsic curvature is defined implicitly via y,^-, oc, j3', and their derivatives through the 
ADM time evolution equation for y,y (see eqn. JA.Sh . 



6.2.2 Painleve-Gullstrand coordinates 

PairJeve-Gullstrand coordinates ICoo02l section 3.3.3] are interesting because they have a flat 
three-metric y,j = 5ij. However, it seems empirically that their time evolution leads to higher 
discretisation errors. The spacetime is described by 



Tij 



I2M 



(6.35) 

(6.36) 
(6.37) 
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(6.38) 

where the parameter M is the mass of the black hole. The auxiliary quantity q' is defined as 
q' = x'^/p.l do not allow for spin, p is the usual radial coordinate with = + + z^. 



6.2.3 Harmonic coordinates 

Harmonic coordinates follow from the coordinate conditions Ox^^^ = ICoo02l section 3.3.2]. A 
black hole in harmonic coordinates is described by 

m = +(^v + v^ + v^^ q'qi (6.39) 
a = \ , (6.40) 

= a^v^q' (6.41) 

where v = 2M/p, and q' = x'/p. p is the usual radial coordinate with p^ = + + z^, and M is 
the mass of the black hole. 1 do not allow for spin. The extrinsic curvature is defined implicitly via 
Yij, a, I3\ and their derivatives through the ADM time evolution equation for y/y (see eqn. tA.St l. 



6.2.4 Coordinate transformations 

In order to have access to a larger class of initial data and analytic solutions, I allow for a generic 
coordinate transformation to be applied to the initial data.^ This transformation is applied to the 
solution's four-metric, making it automatically covariant. The transformation can e.g. be used to 
rotate a black hole to point the spin into any direction, boost the black hole, deform the black hole 
in various more or less useful ways, and change to a moving coordinate system. 

Such a generic coordinate transformation is defined by the usual x^ = Ti^x^ with the trans- 
formation tensor Ty = dx^ /dx^. I also allow for a coordinate translation x^ = x^ + C^. This 
makes for altogether 16 + 4 = 20 free parameters. In order to allow the coordinate change to be 
prescribed conveniently, I decompose the coordinate transformation into 

T = Roto Slow o Dfrm o Cshn o Shear o Bsto Vel (6.42) 

where 

Rot is a rotation of the spatial components. 
Slow is a slowdown, i.e. a scaling of the lapse, 

Dfrm is a deformation of the spatial components, i.e. a rescaUng of the coordinate axes, 

Cshn is a cushion deformation of the spatial components. 

Shear is a shear transformation of the spatial components, 
^Currently I implemented this only for the Kerr-Schild solution. 
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Bst is a Lorentz boost, and 

Vel is a change to a moving coordinate system. 

Of the above, Bst is not really necessary and could be replaced by a combination of the other trans- 
formations, but I keep it for convenience. The individual transformation operators are defined as 



Roty 


= 5'j-a[n.''j (coslncx-l) 


(6.43) 




+ Cl'j sin27ra 


(6.44) 




with = e,^7f rot'Ya 


(6.45) 




and (X = |rot| 




oXOWq 


= slow 




Dfrmj. 


= diag (df rm'^) 


(6.48) 


Cshnj 


= 6) + \eijk\cshn'' 


(6.49) 


Shear^. 


= 5'j + Ciji^ shear*^ 


(6.50) 


Bstj] 


= r 


(6.51) 


BstQ = Bst 


= bst' y 


(6.52) 


Bst'j 


y2 

= 5': H — bst' bst^ 

' Y+1 


(6.53) 




with y = ^ 


(6.54) 


Velj, 


= vel' 


(6.55) 



where the remaining components of these four-tensors are set to 8%. The lower-case three-vectors 
rot', df rm', cshn', shear', bst', vel', and the scalar slow parameterise these transformations with 
19 (instead of the necessary 16) free parameters. Additionally there are the 4 translation parameters 
C^. 

A rotation is described by a three- vector rot'. Its direction is the axis of rotation, and its length 
the angle, where a length of 1 means one full rotation. A cushion transformation changes a square 
into a cushion- or a barrel-shaped object, where the components of the three-vector cslin de- 
scribe the distortion factors along the diagonals. A shear transformation changes a square into a 
rhomboid, where the components of the three-vector sliear*^ describe the shear factor along the 
corresponding coordinate axis. 

Rotations are useful e.g. to point the spin of an object in an arbitrary direction. Deformations and 
shears can be used e.g. to create coordinates where the contravariant and covariant three-vectors 
differ (for testing purposes). A velocity transformation can be used to cancel effects of a boost, e.g. 
to create a stationary boosted black hole (which differs from a non-boosted black hole). Rotations 
and boosts together form the proper Lorentz transforms. 

From the transformed four-metric I then calculate the ADM quantities, i.e. the three-metric, the 
extrinsic curvature, lapse, and shift. The extrinsic curvature is defined via the time derivative of 
the three-metric using equation | |A.5> . The necessary partial derivatives can easily be calculated 
numerically to any given accuracy, much more accurately than to the grid spacing accuracy. 
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When lapse and shift are later, during the time evolution, determined through elliptic equations, 
one needs boundary conditions for them. I often use the boundary values from the initial data as 
Dirichlet boundary conditions. This makes the method used to initially calculate lapse and shift 
actually important. 



6.3 Multiple black hole data 

6.3.1 Brill-Lindquist data 

Brill-Lindquist data ICoo02l section 3.1.2] can contain arbitrarily many black holes. The black 
holes are described by their coordinate location and by a mass parameter. This solution is also 
only valid on a single time slice. It is time-symmetric and conformally flat: 

Yij = M>^^ij (6.57) 

Kij = (6.58) 

oc = 1 (6.59) 

i3' = (6.60) 

where x„ are the positions of the black holes, and ]in are their mass parameters. 

If there is only one black hole, then the convention to write ii/2 instead of 2^ means that the 
event horizon is at r = )i/2 instead of at r = 2fi. In this case, the black hole still has the mass 

M = 

6.3.2 Superposed Kerr-Schild data 

For multiple black hole runs I prefer superposed Kerr-Schild data. They were to my knowledge 
first proposed by Mazner et al. LMHS99II and have recently been refined by Moreno et al. |,MN S02 1 . 
They follow a rather intuitive approach. One views a single black hole as the sum of a flat space 
metric and a black hole metric contribution. This view is motivated by the form of equation <6.23t 
in which the Kerr-Schild four-metric can be written. One then combines two black holes by adding 
two different black hole contributions to a flat space metric. This combined metric does, however, 
not satisfy the constraints any more, so that a constraint solving step has to follow. 

The black holes can be combined in two ways. In the first way, one combines the three-metrics 
and the extrinsic curvatures; in the second way, one combines the four-metrics and their time 
derivatives. Both ways are equally "valid" in principle. Note, however, that the second way breaks 
down when the black holes are close to each other. It can be convenient to attenuate the combined 
solutions as one gets close to one of the black holes, regaining a single black hole solution near the 
individual holes. 

When superposing the three-metrics and the extrinsic curvatures of n black holes, one uses 



Yii = 5, 



•1 



(«) X ' 



(6.61) 



^The combined four-metric then does not have a {— , +, +, +) signature any more. This cannot happen when superposing 
the three-metrics and extrinsic curvatures. 
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13' 



in) 



a 



in) 



{n)i 



(6.62) 
(6.63) 
(6.64) 



The quantities -^"^ are the corresponding quantities from the to-be-superposed black holes. Note 
that the way in which the extrinsic curvatures are superposed means that the superposition of two 
black holes at the same position does not give a single black hole with the combined mass, and 
does not satisfy the constraint equations any more. 

When superposing the four-metric and its time derivative, one uses 



guy 



(n) 



(6.65) 



and then calculates three-metric, extrinsic curvature, lapse, and shift from this four-metric and its 
time derivatives. In both cases, one can also use a different lapse and shift than calculated here, 
but this will lead to a different extrinsic curvature with the second method. When attenuating, one 
attenuates the contributions from the individual black holes. 

As mentioned above, the resulting spacetime will in general not satisfy the constraints. One has 
to use the combined metric as background metric and initial guess to explicitly solve the constraint 
equations (see chapter]^. I usually use the conformal factor \p resulting from the superposition 
and vector potential Vj = as boundary conditions for this. It is also possible to use a Robin 
boundary condition for xp. In order to verify that this procedure does indeed lead to black holes, 
one has to locate the apparent horizons (see section lT^ . 
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The computer code that is used by us numerical physicists is of utmost importance to us. It is 
our experimental setup, the testbed for our ideas, the centre about which our group assembles 
in the morning, and the equipment that we check upon even on weekends. Small codes might 
be written on rainy afternoons or on lazy weekends; they come a dime a dozen and are quickly 
forgotten. The large codes take several people and several months or years to assemble and fine- 
tune, and correspondingly one has to put a lot of effort into keeping a code maintainable and 
understandable to others. 

Yet a code is never, at least not to us numerical physicists, an end in itself. It is a tool; it is an 
important one, but yet only a tool. 



7.1 The Tiger code 

I created and used the Tiger code^ to test the ideas presented in the previous chapters. It descended 
from the Maya code IPenI , written in 2000 at Perm State under Pablo Laguna, which in turn de- 
scended from the Agave code fCol], which in turn was one of the results of the Binary Black Hole 
Coalescence Grand Challenge |A11|. 

The Tiger code is written in the Cactus framework lABD"'"01llABG"'"()TllCacl . which relieves the 
programmer of many merely software-engineering related problems, and which is also supposed 
to make easier the sharing of code modules. At the time of this writing, the sharing still has 
to happen. The main obstacles seem to be rather mundane things such as e.g. differing variable 
names. A promising standardisation effort for ADM-like evolution codes was started in the spring 
of 2002 by the home institution of Cactus, the Albert-Einstein-Institut in Golm near Potsdam. 

When writing a numerical code for a complicated system of equations, there is always the issue 
of how to translate the equations into code. One can either create the code semi-automatically with 
a symbolic algebra package such as Maple fMapf, or one can code the equations by hand. Both 
approaches have their advantages. I chose to code the equations by hand. In doing that, I strove 
for clarity, relying on the compiler to produce efficient code. Appendix IB.4I shows some example 
code. 

The Tiger code has the standard structure found in time evolution codes. (This structure is also 
mandated by Cactus.) It consists of five large components: initial data, constraint solving, time 
evolution, boundary conditions, and analysis routines. Each of these components is described 
below. 



^TGR, pronounced "Tiger": the Tubingen General Relativity Code. (Codes apparently have to have names.) The pronun- 
ciation "Tiger " was suggested by Gabrielle Allen. 
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7.2 Initial data 



In order to facilitate exchanging initial data routines and initial data themselves, it is customary 
in numerical general relativity to create initial data in the ADM variables (see chapter |6), i.e. for 
the three-metric y,^- and the extrinsic curvature While these are sufficient as initial data, it is 
also often customary to additionally provide the initial lapse oc and shift j3'. With lapse and shift 
included, the whole four-metric gf^v can be reconstructed on the initial time slice. 

The Tiger code contains initial data routines for several analytic solutions and background data. 
Among these are the Minkowski spacetime (flat space, see section l6.1.1> , weak Bondi waves (lin- 
ear planar gravitational waves, see section l6.1.21 , the background data for Brill waves (see section 
I6.1.4> , and black holes in Kerr-Schild coordinates (see section l6.2.1> . Painleve-Gullstrand coordi- 
nates (see section l6.2.2t , and harmonic coordinates (see section |6. 2. 3> , as well as Brill-Lrndquist 
black holes (see section l6.3.1> . The Kerr-Schild data can have generic coordinate transformations 
applied (see section l6.2.4> , and they can also be used as background data to superpose black holes 
(see section l6.3.2> . 

The backgrotmd data mentioned above do not satisfy the constraints. They can be used as 
background and as initial guess to solve the constraint equations numerically (see below), which 
is necessary to obtain a solution to Einstein's equations from them. 

Sometimes it is convenient to create initial data e.g. in spherical or cylindrical coordinates, and 
then transform these into Cartesian coordinates. For this, the initial data are still calculated at 
the grid points forming the Cartesian grid, but the tensors have their components in another co- 
ordinate system. A coordinate transformation into Cartesian tensor components is then a linear 
transformation according to the usual T' = {dx' /dxi)TK 

After creating initial data in the ADM variables, these are converted to the TGR variables (see 
chapterUJ, which form the primary variables in the Tiger code. 



7.3 Constraint solvers 

After creating initial data in the TGR variables, and after every time step, these data are used 
as background and as initial guess to enforce the gauge and the constraints, and to calculate the 
lapse and shift. This is the part that distinguishes the Tiger code from other codes that evolve 
unconstrained and without fixing the gauge. 

Enforcing the gauge condition for F,, enforcing the constraints, and calculating lapse and shift 
requires solving nonlinear coupled elliptic equations. This is currently done using the thorn 
TATelliptic in Cactus. This thorn is an interface to generic nonlinear elliptic solvers. The only 
currently available reasonable solvers are TATPETSc, which is an interface to the PETSc library 
iBBG+OlllBGMSgyilBGMgoTl , and TATMG, which is a full approximation storage multigrid solver 
IWes92l that is currently being written by me. 

Enforcing the gauge condition F, yields the traceless conformal metric hjj. In order to calculate 
Yij from this, one has to choose trace y,^ such that dety,^ = 1. This is a nonlinear equation that 
does not involve derivatives. I solve it with the zriddr routine of Numerical Recipes ll-'TVE92l . 
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7.4 Time integration 

The evolution equations are the equations for the time derivatives of the primary variables. In the 
Tiger code, these are z/), y,j, and A,j (see chapter|4j. Evaluating these equations is straightforward 
but tedious (see appendix IA.2.21 eqns. {A.21} . and One usually introduces the 

Christoffel symbols F^^j. and/ or Fyj. and the Ricci tensor R/y as explicit intermediate quantities. 

The time evolution equation for i/) is not strictly necessary during time evolution, but is there 
nonetheless for analysis purposes. The Tiger code contains additionally the time evolution equa- 
tions for the gauge variables K and F,, which are needed to determine the lapse a and the shift /3' 
(see appendix lA.2.2> . 

The time evolution equations for the Hamiltonian and momentum constraints are not imple- 
mented. One could theoretically use them to verify the Bianchi identities at run time, and thus 
gather an additional measure for discretisation errors. 

The time integrator proper uses the iterative Crank-Nicholson scheme ITeuOOl (see appendix 
IB.2> . This is an explicit second-order scheme that can be used to introduce some diffusion in 
order to obtain a stable discretisation for advection terms. As the name suggests, this diffusion is 
introduced via additional iterations. 1 often run with just a single iteration, so that the scheme is 
identical to the midpoint rule (see appendix lB.2t and does not add any artificial diffusion. Instead, 
I add explicit artificial diffusion terms to the time evolution equations (see appendix lB.2.1t . 

7.5 Boundary conditions 

Boundary conditions are necessary for enforcing the constraints and gauge conditions, and for the 
time evolution. 

7.5.1 Outer boundary 

For the constraint, gauge condition, and coordinate condition solvers, there are Dirichlet and 
Robin boundary conditions available on the outer boundary. Dirichlet boundaries are kept con- 
stant while solving, but they do not have to be constant in space. Robin boundary conditions 
enforce a certain falloff towards infinity (see section l5.1.2> . They are readily availabe in Cactus 
|Ca3. 

For the time evolution, there are Dirichlet and radiative boimdary conditions available. Dirich- 
let boundary conditions are kept constant in time, but they do not have to be constant in space. 
Radiative boundary conditions assume that the evolved quantity is an outgoing spherical radial 
wave (see section l5.1.2> . They, too, are availabe in Cactus. 

Additionally, certain symmetry conditions can be enforced at the boundaries, e.g. to restrict the 
simulation domain to a quadrant or an octant of space, or to enforce periodicity. 

For the robust stability tests (see section lS^ . noise can be introduced at the boundary. 

7.5.2 Excision boundary 

For the constraint, gauge condition, and coordinate condition solvers, there are only Dirichlet 
boundary conditions available on the excision boimdary. These Dirichlet boundaries are kept 
constant while solving, but may vary in space. 
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For the time evolution, there are Dirichlet and extrapolation boundary conditions available on 
the excision boundary. Dirichlet boundaries are kept constant in time, but may vary in space. 
The extrapolation boundaries apply an n-th order polynomial extrapolation (n G [0 ... 3]) towards 
the excision boundary. They are applied along the normals of the excision boundary. The high 
order of extrapolation is necessary because it is applied close to the singularity, where the shape 
of the extrapolated functions is not anjrwhere near flat. It is quite possible that a non-polynomial 
extrapolation (e.g. in terms of 1 /r) would perform better. 

7.6 Analysis routines 

After every time step, the ADM variables are calculated from the TGR variables. The ADM vari- 
ables by themselves are important analysis quantities, because they allow comparisons to analytic 
solutions, or to results from other codes. However, because the coordinate system that is used to 
represent the ADM variables is generally not known, the three-metric and the extrinsic curvature 
alone do not provide much insight into the result of a simulation. While they are, together with 
the lapse and shift, sufficient to extract the four-metric and therefore theoretically all information 
about the spacetime, these quantities are basically impossible to interpret when visualised on their 
own. 

The Tiger code provides, as analysis quantities, also all the intermediate quantities that are cal- 
culated during the time evolution. These are the three-Ricci tensor, the time derivatives of the TGR 
variables, the constraints, the gauge variables and their time derivatives, and the gauge violations. 
It also provides the time derivatives of the ADM variables, and the constraints as calculated from 
the ADM variables. Additionally it calculates the four-Weyl tensor, which can be used to extract 
information about gravitational radiation from a simulation, although 1 did not attempt this in my 
test runs. 

Apparent horizons ITho931 chapter 4], ITho96UShi97llSTJ00llHCM0nllSHM0nilSA02l are the only 
locally (in time) detectable structures that are present in vacuum spacetimes. Event horizons are 
global structures and cannot be detected during a time evolution. The Tiger code can locate and 
track multiple apparent horizons. It reports the area and irreducible mass of the horizon. If the 
horizon is isolated IDKSSI . it calculates also the spin and total mass. 
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In this chapter I present results from example runs with the Tiger code, using the formalism 
laid out in the previous chapters. 1 begin with some convergence tests demonstrating that the 
Tiger code is second-order convergent apart from artificial diffusion, and also showing the typical 
magnitude of the discretisation errors. A section about runs with added noise then supports my 
claim that the formalism and its implementation are robustly stable in the sense of Szilagyi et al. 
feGBWOO SSW02J . Two sections on Brill waves and Kerr-Schild black holes show that the code 
works well in the presence of strong field dynamics and black holes. 

8.1 Convergence tests 

8.1.1 Static and stationary tests: black holes 

Here 1 present a series of convergence tests for highly nonlinear but static or stationary spacetimes. 
By calculating the time derivatives of several quantities, essentially the whole right hand side 
evaluation subsystem of the code is tested. The tests here contain no time evolution (but all the 
following tests do). All time derivatives should converge to zero to second order as the numerical 
resolution is increased. 

As test cases 1 use black holes in several coordinates, namely Kerr-Schild (see section (6.2. It , 
Painleve-GuUstrand (see section l6.2.2> , and harmonic coordinates (see section l6.2.3> . In all cases, 
the mass of the black hole is M = 1. For Kerr-Schild coordinates, the spin is = 0.9, and for the 
other two cases, the spin is zero. These are stationary or static solutions, meaning that the analytic 
values of the time derivatives are all zero. The time derivatives of certain primary quantities are 
shown, as well as some constraints. All shown quantities should be zero. The absolute values of 
the numerical values along the x axis are plotted, where the results from the finer resolutions have 
been scaled by factors of 4 and 16, respectively. The dips in the graphs indicate places where the 
numerical error vanishes because it changes sign. 

The fact that the graphs for the different resolutions in figures (0118.21 and l8.3l do overlap indi- 
cates second order convergence. Nevertheless, the errors of some of these quantities are surpris- 
ingly large. Note that the coarsest resolution that 1 present here is dx = 1/16, which is already 
considered to be a rather fine resolution for a typical black hole collision simulation. 

The L2 norms of the convergence factors for the resolutions dx = 1/32 and dx = 1/64 and 
for the three coordinate systems (cs) Kerr-Schild (KS), Painleve-GuUstrand (PG), and harmonic 
coordinates (he) are given below. The fact that the factors are so close to 4 indicates near-perfect 
second order convergence of the code: 



cs 




dtK 


dthxx 


dtAxx 


H 


Mx 


KS 
PG 
he 


3.99554 
3.99958 
3.99976 


4.00827 
4.00089 
3.99453 


4.00078 
3.99852 
(4.59489) 


4.00763 

4.0006 

4.00214 


3.99803 

(n/a) 

3.99408 


3.99002 
3.99987 
3.99923 
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The convergence factor for the Hamiltonian constraint H for Painleve-Gullstrand coordinates 
is ill-defined because the metric is flat, and hence the Hamiltonian constraint is exactly zero ev- 
erywhere. The convergence factor for dthxx for harmonic coordinates is larger than 4 because the 
discretisation error changes sign at x = 2 and hence the convergence factor is ill-defined there. 
Excluding this grid point, the L2 norm of this convergence factor is 4.00058. 
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;ure 8.1: Convergence test with Kerr-Schild coordinates. The graphs show the errors in the time 
derivatives of various quantities for three resolutions. The errors for the finer reso- 
lutions have been scaled by the factors 4 and 16, respectively. The coinciding graphs 
show that there is second order convergence towards zero. 
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Figure 8.2: Convergence test with Painleve-GuUstrand coordinates. The graphs show the errors 
in the time derivatives of various quantities for three resolutions. The errors for the 
finer resolutions have been scaled by the factors 4 and 16, respectively. The coinciding 
graphs show that there is second order convergence towards zero. 
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;ure 8.3: Convergence test with harmonic coordinates. The graphs show the errors in the time 
derivatives of various quantities for three resolutions. The errors for the finer reso- 
lutions have been scaled by the factors 4 and 16, respectively. The coinciding graphs 
show that there is second order convergence towards zero. 
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8.1.2 Dynamic linear test: weak Bondi wave 

As a test of the d5mamic behaviour in the linear regime I use weak Bondi waves (see section l6.1.2t . 
Because a periodicity boundary condition is not possible when solving the elliptic constraint and 
gauge equations, I decided to use standing waves with Dirichlet boundary conditions instead. It is 
also possible to use Dirichlet boundaries for a travelling wave, but they would continuously inject 
the analytic solution into the simulation domain, which I want to avoid. By putting the Dirichlet 
boundaries at nodes of the standing wave, the boundary values are constant in time, which is a 
stronger test of the code. 

I use K = and F, = as gauge conditions. This is consistent with the analytic solution. 

I use an effectively one-dimensional simulation domain, which corresponds to having a trans- 
lational symmetry in two directions. My simulation domain extends only into the z direction with 
z e [—0.5; +0.5]. I use a wave length L = 1 and an amplitude of A = 10^^. I choose a resolution 
of dx — 1/100 and add artificial diffusion with a coefficient Csm = 1/4 (see appendix lB.2.1> . I use 
the midpoint rule (see appendix lB.2^ for time integration. 

Figure l8^ shows the result of this simulation. The phase of the wave stays correct, which is to 
be expected for a standing wave. The amplitude decreases with time, which is also to be expected 
due to numerical and artificial diffusion. The quantities gzz, Fz, and the constraints H and Mz are 
all correct up to floating point accuracy. 

The amplitude decay is exponential in time. It can be empirically described by the expression 
A{t) = A(0) e^^* with a resolution and wave length dependent decay rate t. The table below 
shows decay rates for several example resolutions and wave lengths: 



dx 


L 


A{t = 1)/A{t = 0) 


T 


t/to 


1/100 


1 


0.906 


0.0987 (=: To) 


1 


1/200 


1 


0.952 


0.0493 


0.500 


1/100 


1/2 


0.673 


0.3964 


4.017 


1/100 


1/4 


0.207 


1.5764 


15.98 



This table indicates that the decay rate t is proportional to the resolution dx and inversely pro- 
portional to the square of the wave length L. 

This rate can be explained by the kind of artifical diffusion that I use: Given initial data with 
A <C 1 as described in section l6.1.2l gauge conditions as described above, and including artificial 
viscosity with Csm 7^ as described in section lB.2.11 the set of equations governing the time 
evolution of 

dtgxx = -2Kxx + CsMdxdzzgxx (8.1) 
dtKxx = -^dzzgxx + CsMdxdzzKxx + 0{A^) . (8.2) 

Neglecting the terms in O(A^), these equations can be combined into 

dttgxx = dzzgxx + 2CsMdxdtzzgxx + 0{dx^) . (8.3) 
Neglecting again the terms in 0{dx^), the ansatz 

gxx{t,z) = exp{—k^ Csudx t) sin{kt) cos{kz) , (8.4) 
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Figure 8.4: Performance test with a weak Bondi wave for the resolutions dx = 1/100 (top left) 
and dx = 1/200 (top right). Shown are the wave forms at certain times during the 
evolution. The bottom graph compares the amplitude decay for the two resolutions. 
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which is motivated by the form of the initial data, solves this equation. That means that the decay 
rate t caused by the artificial diffusion should be 

CsM dx 

T = — jj — (8.5) 

with L = 1/k, and this is indeed the case in the table above. It is therefore reasonable to assume 
that the amplitude decay is due to the artificial diffusion. 

Although an amplitude loss of 10% per crossing time seems large, it is actually acceptable for 
a binary black hole collision simulation. Scaling the resolution to typical values for such a run, I 
arrive at a simulation domain with x' G [0;20], a resolution of dx = 1/5, and a wave length of 
L = 20. This wave length is close to the quasi-normal mode of a black hole with a mass of M = 2, 
the result of a collision of two M = 1 black holes. In this case, it takes a time T = 20 for the 
numerical amplitude to decrease to 90% of its real, physical value. This is about the time scale 
in which the wave reaches the outer boundary. That means that the code is not perfect, but the 
performance^ is acceptable for a realistic run. 



use the term "performance" as measuring the quality of the numerical result. This performance depends on the mag 
nitude of the numerical errors, not on the speed of an implementation on some hardware. 
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8.1.3 Dynamic nonlinear test: gauge pulse 

As a test of the dynamic behaviour in the nonlinear regime I present a nonlinear gauge pulse (see 
section l6.1.3t . For the same reason as in the previous section, a periodicity boundary condition is 
not possible. Again, 1 therefore use a standing gauge pulse with Dirichlet boundary conditions 
that are constant in time. 

1 use the analytically known (time and space dependent) values of K and f , as gauge condi- 
tions. This is permissible, because these are gauge quantities only, and one is free to specify the 
gauge condition in advance — indeed, one has to. It is certainly possible to use a different gauge 
condition, but one then cannot compare the result to this analytic solution any more. 

I use an effectively one-dimensional simulation domain, which corresponds to a translational 
symmetry in two directions. My simulation domain extends only into the z direction with z E 
[—0.5; +0.5]. 1 use a wave length L = 1 and an amplitude of A = 1. I choose a resolution of 
dx = 1/100 and add artifical diffusion with a coefficient Csm = 1/4 (see appendix IB. 2. 1> . 1 use 
the midpoint rule (see appendix IB. 2> for time integration. Except for the wave amplitude, these 
parameters are the same as in the previous section. 

Figure (831 shows the result of this simulation. The phase stays correct, which is to be expected 
for a standing wave. The amplitude does not decrease, it stays constant in time. The simulation 
result does not deviate from the analytic solution in any significant manner. 

This test case contains a pure gauge wave, and is therefore a pathologically well suited^ case for 
my formulation. Because I explicitly enforce the gauge condition at every time step, the nontrivial 
part of the evolution for these data is prescribed. The gravitational wave degrees of freedom are 
not prescribed, but those stay zero all the time. 



^The term "pathologically well suited" was coined by Ed Seidel in spring 2002 in Mexico. 
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Figure 8.5: Performance test with a gauge pulse. The top graph compares the initial wave form and 
the wave form after 100 crossing times. The bottom graphs show the wave amplitude 
(left) and constraint and gauge violations (right) vs. time. 
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Figure 8.6: Robust stability test. Plotted are the violations of the constraints and of the gauge con- 
dition versus time for the resolutions dx = 1/4 (left) and dx = 1/8 (right). 



8.2 Stability test 

Szilagyi et al. ISGBW00IISSW02I define the notion of robust stability of a code. This definition is not 
an analytic, but rather an experimental one, which means that it can be tested rather easily. The 
basic idea is to add noise to the initial data, and also to add noise to the values provided by the 
boundary conditions. A code that behaves gracefully, i.e. does not show exponential growth, is 
called robustly stable. 

The robust stability test, as proposed by Szilagyi et al., comes in four stages of increasing diffi- 
culty. Stages I and II include periodicity, which is not a valid boundary condition for my elliptic 
equations. I therefore concentrate on stage III, which uses a rectangular domain without periodic- 
ity. Stage IV requires a spherical outer boundary, which I did not test. 

For this test, I use Minkowski flat space (see section l6.1.1> as background to which the noise is 
added. I use K = and F, = as gauge conditions and Dirichlet boundary conditions. 

I run this test in two configurations with different resolutions dx = 1/4 and dx = 1/8. I use a 
cubic box with a length of L = 4, and run the code up to T = 400, which is equivalent toT/L = 100 
crossing times. I use a noise amplitude of A = 0.1 on the initial data and on the boundaries. I add 
an artificial diffusion with a coefficient of CgM = 0.1 (see appendix lB.2.1> . I use the midpoint rule 
(see appendix lB.2> for time integration. 

The result of this test is that the formulation of Einstein's equations as implemented in the Tiger 
code is robustly stable.'' The constraints and the gauge violations, which are shown in figure 1531 
do not increase with time after an initial transient. Similarly, all components of the three-metric 
and of the extrinsic curvature stay bounded. 

I assume that the artificial diffusion is the reason for the fact that there seems to be no growth 
at all. I assume that the levels at which the individual metric and extrinsic curvature components 
remain are defined by the noise on the boimdary, which inserts energy into the system, and the 
artificial viscosity, which removes energy. 

Note that the finer run is not just a higher-resolution version of the coarser run. The two space- 
times are different, in that the finer run contains noise with a higher spatial frequency. For that 

longer run time would have been preferable, but was not possible because of the large computational resource re- 
quirements. 
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reason the errors in the fine run cannot be expected to be smaller than those in the coarse run. The 
two runs do not form a convergence test. 
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8.3 Brill wave 



A Brill wave (see section 16. 1.4> tests the dynamic nonlinear behaviour in a realistic situation. Al- 
though the initial data that 1 use are axially symmetric, I perform a full three-dimensional evolu- 
tion. Due to the Cartesian grid in the simulation domain the setup loses its axial symmetry. 

There is no complete analytic solution for the Brill wave initial data; the conformal factor has 
to be determined numerically from the Hamiltonian constraint. In addition, there is no analytic 
solution at all for the evolution of Brill waves at times t ^ 0. It is thus difficult to gauge the 
correctness and performance of a code. 

The Albert-Einstein-lnstitut in Golm also has a code, called Einstein code, that solves Einstein's 
equations. This code is publicly available from the Cactus web pages ICacl . 1 used both their and 
the Tiger code, and tried to compare the results for Brill waves. This is unfortunately very difficult, 
because the codes use different gauge conditions during the evolution. While both codes can use 
maximal slicing, i.e. K = 0, they differ in their shift conditions. The Tiger code uses a metric gauge 
condition F, = const or F, = 0, and derives a shift condition from that. The Einstein code cannot 
impose a gauge condition on the metric, and one has to specify the shift condition directly. For 
Brill waves one usually uses normal coordinates, i.e. /3' = 0. This difference leads to very different 
time evolutions, although the spacetimes should (barring implementation errors) be identical up 
to the discretisation error. 

I simulate Brill waves with initial data as described in section l6.1.4l 1 use an amplitude A = 4, 
which makes the Brill wave highly nonlinear, but still subcritical, i.e. it does not form a black hole. 
1 also use C = and w = 1. I use maximal slicing, i.e. the gauge condition K = on the extrinsic 
curvature. 1 use an iterative Crank-Nicholson time integrator with 2 iterations after the initial 
Euler step (see appendix lB.2^ . 

1 run the Tiger code with two different gauge conditions. The first, which I call here "F, = const", 
keeps the gauge variable F, constant in time at the values from the initial data. The second, called 
"Fj = 0", sets F, to zero at all times. Thus only the case F, = const uses the same initial data as 
the Einstein code. The case F, = evolves a different (but closely related) spacetime. This gauge 
condition has the advantage that it, together with K = 0, enforces a Minkowski metric when the 
spacetime is flat (and with suitable boundary conditions). 

It turns out that the cases F, = const and F, = show very similar behaviour. This indicates that 
my way of enforcing the gauge leaves the physical degrees of freedom mostly unchanged. This 
is also confirmed by the time evolution of the ADM mass, as shown in figure IS^I Strangely, the 
behaviour of the Einstein code using a zero shift differs greatly, as can also be seen in the following. 

Traditionally, one of the most interesting quantities to look at in maximal slicing is the lapse a, or 
its minimum. With maximal slicing, the lapse drops to zero near a singularity; thus the minimum 
of the lapse indicates heuristically whether a singularity is forming. The Brill wave that I look at 
is subcritical, i.e. there is no singularity in the spacetime. Consequently, although the lapse drops 
initially, it later "recovers" and ends up as a: = 1 everywhere at late times. 

These time evolutions of the minimum of the lapse are shown in figure 15^ for four resolutions 
from dx = 1/4 to (ix = 1 /8, run with the Tiger code. Finer resolutions would have been desirable, 
but would have required much longer run times. (The PETSc elliptic solver (see section I7!5l that 
1 used does not scale well on multiple processors.) As the resolutions are, they form only a weak 
convergence test. The same figure also shows a run with the Einstein code with a resolution of 
dx = 1 /6. Figures l8^ and l8 . 1 01 show the time evolution of the lapse along the x axis. 

All three evolutions show qualitatively the same features: the lapse first drops sharply down 



61 



8 Results 



ADM mass ADM mass 




2 4 6 8 10 5 10 15 20 

t t 

Figure 8.7: Brill wave runs: ADM mass vs. time for different gauge and coordinate conditions. 
Note the different time scales. 



from a = 1, bounces back up two times, and then drifts back towards a = 1. At late times, the Brill 
wave has radiated away, and the spacetime becomes flat. The largest and maybe most puzzling 
difference is that it takes about twice as much coordinate time until the lapse has recovered for the 
zero shift coordinate condition than for the two other gauge conditions. 

This difference in coordinate time seems to be caused by a combination of two effects. First, 
due to the different shift values, the lapse near the origin has different values. This is because the 
shift enters into the dfK equation <A.23> that determines the lapse. (The shifts are shown in figure 
18.111 ) A smaller lapse leads to a larger ratio between coordinate time and proper time. However, 
this explains only a part of the difference, as can be seen when integrating proper time along the 
geodesic formed by the origin (see fieure l8J2l . 

The second reason is that the lapse is a pure coordinate quantity, and by comparing the lapse 
one cannot directly make a statement about the location of gravitational waves. That is, the Brill 
wave might have radiated away before the lapse has fully recovered. A better (but still not perfect) 
quantity to look at is gi/y along the x axis. This is a metric component which is orthogonal to the 
direction in which the Brill wave radiates, and hence is connected to the transverse gravitational 
wave degrees of freedom in the system. 

Figures f8. 131 and 18. 141 show the evolution of gi/y vs. time. One can see that for the gauges F; = 
const and F, = there are two wave trains that leave the simulation domain at x = 6 at about 
t = 10. With the coordinate condition 13' = 0, there are also two wave trains, but it is clear from 
the graph that they arrive at x = 6 substantially earlier than att = 20. Hence there was not really 
a factor of two difference between the coordinate times at which the Brill wave has radiated away 
to begin with. 

Figure ISTSI compares then initial (at t = 0) and late time (at t = 10 and f = 20, resp.) behaviour 
of several quantities, namely the conformal factor xp, the Ricci scalar R, and the transverse metric 
component gyy. All quantities are plotted vs. the x coordinate. 

By construction, the conformal factor ip is initially the same for F, = const and [3' = 0. With 
K = and ji' = 0, it is Dfi/) = (see eqn. l lA.21t '), so that the final xp is the same as the initial one in 
this case. For K = and F, = 0, flat space has a Minkowski metric, so that xp = lin this case. (This 
state has not yet completely been reached att = 10.) Note that the cases F, = const and ^' = go 
to different static coordinate systems of flat space. 
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8.3 Brill wave 



F;=const FpO 




5 10 15 20 

t 



Figure 8.8: Brill wave runs: Minimum of the lapse a vs. time for different gauge and coordinate 
conditions. (This is always the value of the lapse at the origin.) All cases use maximal 
slicing, i.e. K = everywhere. Note the different time scales. 
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8 Results 



Evolution of the lapse for Fi=const 




Evolution of tfie lapse for Fi=0 




Figure 8.9: Brill wave runs: Lapse a vs. radius (along the x axis) and time for different gauge 
conditions. Note that also K = everywhere. 



64 



8.3 Brill wave 



Evolution of the lapse for beta'=0 




Figure 8.10: Brill wave runs: Lapse oc vs. radius (along the x axis) and time for zero shift. Note that 
also K = everywhere. 
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8 Results 



Evolution of the shift for Fi=const 



beta 
0.1 





Evolution of the shift for Fj=0 




Figure 8.11: Brill wave rims: Shift component 13-^ vs. radius (along the x axis) and time for different 
gauge conditions. Note that also K = everywhere. 
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8.3 Brill wave 



Proper time at the origin Proper time at the origin 




Figure 8.12: Brill wave runs: Proper time vs. coordinate time at the origin for different gauge and 
coordinate conditions. 

The Ricci scalar should be zero initially, due to the Hamiltonian constraint and time symmetry. 
The Ricci scalar is nonzero because of discretisation errors. The final Ricci scalar should be zero 
for the Fj = case. 

Again, by construction the initial values of gyi/ are the same for f,- = const and 13' = 0. The case 
F,- = differs, because it has a different gauge condition enforced onto it. At late times, gyy tends 
to its Minkowski value of one for F; = 0. 

A "real" comparison between the two F/ = const and 13' = runs would involve finding a 
four-coordinate transformation between the two four-metrics. The two spacetimes are identical 
if and only if such a coordinate transformation exists (up to nimierical errors). This would be a 
difficult undertaking, not only because it involves a huge amount of data (namely two complete 
spacetimes). 1 know of no proven numerical method to find such a transformation. The Lazarus 
project IBCL02I tries to solve a similar problem; they compare e.g. curvature invariants. 
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8 Results 



1 

0-9 

Syy 




1 

0-9 

Syy 




Figure 8.13: Brill wave runs: Metric component gyy vs. radius and time for different gauge condi- 
tions. Note that also K = everywhere. 
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8.3 Brill wave 



Evolution of the metric component g„ for beta'=0 




Figure 8.14: Brill wave runs: Metric component gyy vs. radius and time for zero shift. Note that 
also K = everywhere. 
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8 Results 
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Figure 8.15: Brill wave runs: Conformal factor xp, Ricci scalar R, and metric component gyy vs. 

radius along the x axis att = and at f = late. Note that for all shown cases a = 1 and 
13' = 0. The initial Ricci scalar R should be identically zero; it shows the magnitude of 
the discretisation errors. 
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8.4 Kerr-Schild black hole 



8.4 Kerr-Schild black hole 

One of the most basic tests of a nonlinear evolution code is a static or stationary black hole. I 
already tested the right hand side of the evolution equations in various coordinate systems in sec- 
tion lS.l.Il above. Here I want to continue these tests by performing time evolutions with different 
numerical configurations, i.e. with different resolutions, outer boundary locations, and different 
enforced grid symmetries. 

I consider the following configurations: 



title 


resolution 


spin 


location of 
outer bnd. 


symmetries 


std 


1/4 





4 


octant 


dx8 


1/8 





4 


octant 


ob8 


1/4 





8 


octant 


full 


1/4 





4 


none 


spin 


1/4 


1/2 


4 


octant 



In all cases, the mass of the black hole is M = 1 . I excise a region with radius r = 1 around the 
singularity. I use Dirichlet inner and outer boimdary conditions, and I add artificial diffusion with 
a coefficient of Csm = 0.1 (see appendix lB.2.1> . I use an iterative Crank-Nicholson time integrator 
with 2 iterations after the initial Euler step (see appendix IB.2^ . I rim these configuration up to 
t = 100. 

With octant symmetry, only one octant of the spacetime is simulated. The boundary conditions 
at the zero coordinate planes are given by the symmetry conditions. Using these symmetries re- 
duces the memory and run time requirements by a factor of eight, which is considerable. However, 
using or not using these symmetries can influence the stability of the code, and therefore has to be 
tested. 

I also use other, more physical boundary conditions. For these runs, I use radiative outer bound- 
ary conditions, or (third-order, or cubic) extrapolated inner boimdary conditions, or both.^ Other- 
wise, these rims are identical to the std run above: 



title 


outer be 


excision be 


std 


Dirichlet 


Dirichlet 


rad 


radiative 


Dirichlet 


extrap 


Dirichlet 


extrapolated 


extrap-rad 


radiative 


extrapolated 



I show in figure lSn^ the L2 norms of the Hamiltonian and the momentum constraints and of the 
gauge violation. The time evolutions of the above configurations std, dx8, ob8, full, and spin show 
an initial transient, and then settle down close to the initial data, i.e. close to the analytic solution. 
They are all stable in the sense that no growth is visible during the examined time scale. Because 
the code is robustly stable, I do not expect instabilities at later times, although I have no proof for 
that. 

The configurations rad and extrap-rad, which have a radiative outer boundary condition, pick 
up a slowly decaying global breathing mode. This error has an amplitude of about 1% in the ADM 
and apparent horizon masses, as shown in figure This mode has not yet decayed at f = 100. 

^However, the inner boundary condition for the conformal factor xf) was always Diriclilet. 
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8 Results 
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Figure 8.16: Kerr-Schild black hole time evolution: L2-norms of the constraint and gauge 
violations for various configurations and boundary conditions. std=standard, 
dx8=higher resolution, ob8=larger domain, full=without symmetries, spm=with spin, 
rari=radiative outer boundary, exfrap=extrapolated inner boundary 



72 



8.4 Kerr-Schild black hole 
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Figure 8.17: Kerr-Schild black hole time evolution: ADM masses and apparent horizon masses. 

std=standard, dx8=higher resolution, ofc8=larger domain, /w//=without symmetries, 
spm=with spin, rad=radiative outer boundary, extrap=extrapolated inner boundary 
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8 Results 



When comparing simple (std) Kerr-Schild runs with the resolutions dx = 1/4, dx = 1/6, and 
dx = 1/8 at t = 100 along the x axis, I get following convergence factors (cf): 



X 


1.5 


2.0 


2.5 


3.0 


3.5 


cf 


2.26687 


2.71461 


2.74507 


2.73848 


2.73511 



A factor of (1/4^ - 1 /6^)/(l/6^ - 1 /8^) w 2.85714 in the above would indicate second order 
convergence. The ntmibers suggest that the TGR system is able to evolve static black holes in Kerr- 
Schild coordinates with excision in a stable and convergent manner. The stability of the evolution 
is independent of the resolution, the location of the outer boundaries, the symmetries that are 
imposed, and the kind of boundary conditions that are used. 
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9 Conclusion 



In this thesis, I set out to find a way to arrive at some more understanding of the instabilities that 
one encounters in today's black hole evolution simulations. I looked at gauge conditions and clas- 
sified them into gauge evolution and gauge fixing conditions, assuming that gauge modes do play 
an important role in the instabilities. I proposed a gauge fixing condition and an implementation 
for it, and examined it in various test problems and some applications. 

Gauge fixing seems to lead to very stable evolutions. The Tiger code is robustly stable, and has 
also sufficient accuracy to treat more complex configurations. It remains to be shown that this 
system also works for two coalescing black holes. The computational resources required by the 
current implementation do unfortunately not allow such a test; further algorithmic improvements 
(the implementation of a parallel multigrid solver) will probably be necessary for that. 
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9 Conclusion 
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A Equations 



Many of the equations in this appendix are only a repetition of well-known literature. They serve 
mostly as a quick reference, and to define the sign conventions. For other equations, I give a short 
derivation. I restrict myself to the vacuum case everywhere. 



A.l The ADM formalism 

The ADM formalism is one way to transform the Einstein equations into a system of time evolution 
equations. 



A. 1.1 Variables 



The ADM variables IADM62I are the three-metric 7,y, the extrinsic curvature Kij, the lapse a and 
the shift 13'. The four-metric can be written as 







Pi 




Pi 


Tij 







J 

where /3, = YijP^ > arid 0^ = /3y/3A This metric is equivalent to the line element 

ds^ = -oc^df + Yij [dx' + 13' dt^ (^dx' + 0dt^ 

= {-(x^ + Tijl3'l3'^ dt^ + 2jijj3idx'dt + njdx'dx' . 

The inverse metric is then 





/ -1 




1 






P' 


a^r'i - 13^ 






J 



(A.l) 



(A.2) 
(A.3) 



(A.4) 



The three-metric and extrinsic curvature on a time slice are enough to completely specify the 
physics in the whole spacetime. Lapse a and shift j3' can be chosen freely during time evolu- 
tion and determine only the coordinate system. The extrinsic curvature is defined via the time 
derivative of the three-metric by eqn. <A.5> below. 
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A Equations 



A. 1.2 Time evolution 

The time evolution of the three-metric (which is also definition of the extrinsic curvature Kij) is 
given by 

dtYij = -2aKij + rijl3'- + rul3lj + l3'rij,i (A.5) 
The time evolution of the extrinsic curvature is given by 

dtKij = -ViVja + alRij-lKii^+KKij] (A.6) 
+ Kijl5l, + Knl3lj + l5'K,jj 

A. 1.3 Ricci tensor 

The connection coefficients to the three-metric are given by 

rij,k + rik,j - Tjkj (A.7) 

The three-Ricci tensor is e.g. given by 



R.j = r'j,k-4,j + r!fk-^j (A.8) 

It is a bad idea to calculate it in this way numerically, as this would involve taking derivatives of 
derivatives. Instead, one has to explicitely take second derivatives of the three-metric. 

A. 1.4 Constraints 

The Hamiltonian constraint is usually written as 

H = R + K^-KijK'' (A.9) 
The momentum constraint is usually written as 

Mi = \/i{Kij-r,jK) (A.10) 

As these constraints have to be zero for a physical spacetime, one can multiply them by an 
arbitrary nonzero function to arrive at other formulations for the constraints. This especially leaves 
scaling freedoms, so that a statement saying the constraints are smaller than a certain value is only 
meaningful if this is with respect to a certain given length scale. 

One can introduce normalised constraints e.g. as 

= , , , , (A.ll) 

\R\ + \K^ + \KijK'l\ 

m(~) = , . ^\ - (A.12) 

where the absolute value of a vector x' is calculated as | x' p = jijX^xi. These normalised constraints 
have a range of [0; 1] and are thus scaling invariant. They have the disadvantage of being ill- 
defined e.g. in the Minkowski spacetime, and are meaningless in perturbations thereof. 
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A.2 The TGR system 



A. 2 The TGR system 

The TGR system is a conformal traceless ADM (ctADM) formalism that is derived from the ADM 
formalism by introducing two additional scalar quantities xp and K, which represent the determi- 
nant of the three-metric and the trace of the extrinsic curvature, respectively This also leads to two 
new constraints, namely just these conditions. 

A. 2.1 Variables 



The variant of conformal traceless ADM formalism (compare ISN95I |BS99|) that I use has the 
following five variables: 
The conformal factor 



xp' = (detr/;)^/' (A.13) 



The conformal metric 



fij = (A.14) 

where thus dety,-^- — 1 (A. 15) 

The trace of the extrinsic curvature: 

K = Ki (A.16) 
The traceless extrinsic curvature (this is only an intermediate quantity): 

Aj = K,j-^r,jK (A.17) 

where thus A] = (A.18) 

The conformal traceless extrinsic curvature: 

Aj = 4>^Aij (A.19) 

where thus A| = (A.20) 

Quantities with a tilde ~ have their indices raised and lowered using the conformal metric y,y. 
Equations llA.lSt and <A.20> form two additional, algebraic constraints. 



A. 2. 2 Time evolution 

The time evolution equations for the ctADM quantities follow in a straightforward way from their 
definitions, and from the time evolution equations of the ADM quantities. 
The time evolution of the conformal factor is given by 

dtxk = -laxpK+l3'xPj + lxl>l3\ (A.21) 
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A Equations 

The time evolution of the conf ormal metric is given by 



(A.22) 



The time evolution of the trace of the extrinsic curvature is given by 



dtK 



-Ac 



(A.23) 



Here A denotes the covariant Laplace operator V'V;. 

The time evolution of the traceless extrinsic curvature (an intermediate quantity) is given by 



dtK.j - -dt {TijK) 



-V,Vya + (X 



Rij - IKiiK) + KK, 



^'1 



(A.24) 
(A.25) 



1 - 

3 -~ 
1 



2aKij + rijl5[i + rnlilj + pWij,i 



K 



(V,Vya)" + a(R,/^ (A.26) 
-2Ai,A) 
+ Aijl3\^ + Aul3]j + l3'Aij^i 

Here (X^y)^^ denotes the tracefree part of X;y with respect to the physical metric, i.e. Xjy — 

sTzvy'^x,,,, 

The time evolution of the conformal traceless extrinsic curvature is given by 



2xl> --ocxl>K + I3'xl>i + -xpli'i A 



(A.27) 

(A.28) 



-lAuA] 



AhK 



3 'J 



+ xP' [A,jl3'. + Aal3[j + P'Aijj 



-xP^ {ViVja) +aip^ {R 
- laip-^AiiA'j 

+ Aijl^li + Aniilj + (i^Aij^i + ^AijUl, 



(A.29) 
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A.2 The TGR system 

A. 2. 3 Constraints 

The conformal Hamiltonian constraint can be written as 

H = Axp-^ipR-^xP^K^ + ^xp-'AijA'i (A.30) 

The conformal momentum constraint can be written as 

Mi = V^A,y - |i/»^V;K (A.31) 

Note that that these two constraints contain different scahng factors than the ADM constraints 
iA.91 and <A.10> . These two formulations of the constraints cannot readily be compared to each 
other. 



A. 2. 4 Enforcing the constraints 

In order to enforce the constraints, one can extend the system of ctADM variables with a vector 
potential V, for the conformal traceless extrinsic curvature lCoo02i section 2.2.1]. This leads to the 
modified conformal traceless extrinsic curvature 

Ay = A*.+ (Ly),^. (A.32) 

where the operator L is the longitudinal derivative with respect to the conformal metric, defined 
by 

(LV),-y = v,y^ + V;V,-^f,^v'^yt (A.33) 

which makes the gradient (LV);y symmetric and tracefree. 

Given a conformal metric y,j, a mean curvature K and a background conformal traceless ex- 
trinsic curvature A* , one can enforce the constraints by solving the constraint equations for the 
conformal factor t/^ and the vector potential V,-. This vector potential and the background confor- 
mal traceless extrinsic curvature then define the real conformal traceless extrinsic curvature A,y 
through <A.32> . 

A. 2. 5 Gauge condition 

I introduce the traceless conformal metric 

% = yij-\kj^^^yn ■ (A.34) 

The bar " is used to indicate that the indices are to be raised and lowered using the coordinate 
metric b\y 

The gauge variable is 

F, = Wj,^ . (A.35) 
which is close to what Shibata and Nakamura use in ISN95I . 
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A Equations 

The gauge condition should constrain the conformal metric, and the conformal metric already 
has to satisfy the constraint dety,y = 1. Enforcing the gauge condition will change the conformal 
metric, but one has to make sure that it does not change its determinant. By making the gauge 
condition independent of the trace of the conformal metric, one can later change the trace to adjust 
its determinant to one. 



A. 2. 6 Enforcing the gauge condition 

Enforcing the gauge condition on the trace of the extrinsic curvature K is trivial. 

The method that I use to enforce the gauge on the conformal metric is similar to the method 
used to enforce the constraints. I extend the system of ctADM variables with a vector potential Wj 
for the traceless conformal metric. This leads to the modified traceless conformal metric 

hij = h*j + {LW),j (A.36) 

where the operator L is the longitudinal derivative with respect to the coordinate metric, i.e. using 
partial derivatives. It is defined by 

(LW),-^. = Wj, + Wi^j-^5ijWk,k (A.37) 

which makes the gradient (LW),^- symmetric and tracefree. 
The gauge condition is then 

Fi = h*j^j + W,jj+^Wj,j (A.38) 

Given a background traceless conformal metric h*j, one can enforce the gauge condition on the 
metric by solving the gauge equation JA.38t for the vector potential W,. This vector potential and 
the background traceless conformal metric, together with the constraint det fjj = 1 from which 
the trace of the conformal metric is calculated, then define the real conformal metric y,y through 



A. 2. 7 Determining lapse and shift 

The gauge condition implicitly determines the lapse and shift. 

One differentiates the gauge conditions with respect to time, and then uses the time evolution 
equations for the trace of the extrinsic curvature dtK and for the conformal metric dtYij to arrive at 
a set of coupled elliptic equations for the lapse and shift. 

The equation for the lapse a is immediately given by iA.23> . The equation for the shift j3' is lar- 
gish, but not very enlightening, and I omit most of it here. The fact that it is so complicated seems 
to point to the fact that there should be a more elegant gauge condition. The terms containing 
second derivatives of the shift are 

dtF, ^ (A.39) 

= Yik^%+\Yjkli)j-\kjP'^,jk (A.40) 

These terms look sufficiently close to being elliptic that one can hope that the equation is well- 
posed. The niraierical experiment shows that this seems to be the case. 
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B Numerics 



B.l Spatial discretisation 

I discretise the quantities on a Cartesian grid. I approximate the partial derivatives by second 
order centred derivatives. This is among the most simple possible discretisation methods. 



B.2 Time integration 

I integrate in time with an explicit second order integrator, using the iterative Crank-Nicholson 
scheme |TeuOO|. For that I usually use two iterations after the initial Euler step. For some exper- 
iments I use only one iteration after the Euler step, which turns this scheme into the midpoint 
rule, or a second order Runge-Kutta. This is one of the most simple possible second-order time 
integrators. 

Given the ordinary differential equation 



dt 



m 



u[t,f{t)] 



(B.l) 



the midpoint rule uses the time stepping scheme 



(1) 



(2) 



fn-l+h U [tn-l,fn-l] 

/„_i + hu 



1 1 

^{tn-l + tn),-{fn-l+fi'^) 



fn 



f'n 



(2) 



(B.2) 
(B.3) 

(B.4) 



to advance from/„_i to /„ with the step size h. fji^^ is an intermediate quantity. 
Similarly, the iterative Crank-Nicholson method uses the time stepping scheme 



(1) 



(0 



fn-l+hu [tn-l,f„-l 
fn-l +hu 



tn-l + tn),^ {fn-1 +ft^^) 



fn 



lim 



(0 



(B.5) 
(B.6) 

(B.7) 



where the limes is usually only taken up to i = 3. The first step is a simple Euler step. All other 
right hand sides are evaluated at the half-time -|- 
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B Numerics 



B.2.1 Artificial viscosity 

The advection terms in the time derivatives, i.e. those terms involving the shift [3', are unstable 
when discretised as above and integrated with the midpoint rule. I therefore add artificial diffu- 
sion to the right hand sides of the time evolution equations. 1 implement two kinds of artificial 
diffusion, of which one depends on the shift: 

(5t/)AD = L^Cadv^ + Csm)^^'/,. (B.8) 

with < Cadv ^ 1 and < Csm ^ 1- The quantity dx' is the grid spacing. These two viscosities 
are explained below. 

The first kind, ADV, which depends on the shift, transforms the centred differencing advection 
terms into first order upstream differencing terms when Cadv = 1- In one dimension this follows 
as follows: 

/Sar/ ^ Q{^)P ^^''^'\~ ^^-^^ + e(-/3)i3 ^^^^ Y^^"'^^ (B.9) 
= (i3 + |/3|)^ii±^|^ (B.lO) 
-(-^ + 1-/31) ^^-^-^-"^ (B.ll) 

= i3 8.,/+l^ft3,,,/ (B.13) 

where h is the grid spacing. The first line is the first order upstream difference. 1 use the Heaviside- 
function 0(x) which has the value for x < and 1 for x > 0. It is also 2|x|0(x) = x + |x|. 

The second kind of artificial diffusion, SM, is a generic smoothing operator, and is independent 
of the shift. It is necessary is some simulations. 1 found that, without this diffusion, the system 
of equations is unstable against high frequency noise, and the amplification rate depends on the 
grid resolution. This points to either a bad discretisation scheme, or an ill-posed problem already 
at the physical level. 1 did not want to examine this kind of instability any further, so 1 decided to 
remove it by using some suitable artificial diffusion. 

Both kinds of artificial diffusion terms remove high frequency modes, making the simulation 
stable against such noise. They also dissipate energy, dampening the amplitude of gravitational 
waves. They introduce a second order error into the system, turning it into a first order correct 
system only. 



B.3 Elliptic integration 

In order to transform the constraints and the gauge conditions into elliptic equations, 1 add vector 
potentials to the set of equations. That means that enforcing the constraints requires going back 
and forth between different representations, and that introduces discretisation errors. 
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B.3 Elliptic integration 



B.3.1 Variable transformations 

The step of introducing the vector potentials does not introduce errors, as the vector potentials are 
set to zero initially. However, going back to the original system then does introduce discretisation 
errors. This can be demonstrated by using the prototype of a constraint equation 

F, = djhij . (B.14) 

I first introduce a vector potential W„ leading to 

hij = h*j + djW, . (B.15) 

This turns the constraint equation into 

F/ = djh*j + djjW, , (B.16) 

which is now elliptic rn Wj. It can be solved to arbitrary accuracy. Note that F, and F- have 
fundamentally different properties. 

After solving F- for W„ one calculates the new hjj from its definition <B.15> . This requires taking a 
first derivative of W, . In order to evaluate the original constraint F,, one has to take a first derivative 
of hjj, which contains then two successive first derivatives of W,. This is, on the discrete level, 
quite different from the constraint F^', which involves taking one second derivative of Wj. Both 
these steps introduce different discretisation errors. That means, even after solving F^ to very high 
accuracy, F, will still contain an error of the order of the discretisation error. Furthermore, taking 
two successive first derivatives tends to introduce a high frequency noise. 

This problem might be eliminated by using staggered grids. On a staggered grid, two consecu- 
tive first derivatives are the same as one second derivative. 



Double first derivatives 

Using two consecutive first derivatives for solving the equations would be unstable. A second 
order centred first derivative has the stencil weights 

-i +i (B.17) 

while a second derivative has the weights 

+1 -2 +1 (B.18) 

Taking two consecutive first derivatives leads to an effective stencil with the weights 

+ 1 +1 (B.19) 

Stencil JB.19> contains two zeros. If one calculates the derivative of e.g. an even-numbered grid 
point, then only other even-numbered grid points contribute. This decouples odd- and even- 
numbered grid points. If one uses this stencil while solving an elliptic equation, then one obtains 
effectively two independent solutions, namely one for the odd, and one for the even grid points. 
These solutions are coupled through their initial data and boundary conditions only and will differ 
slightly. Going from one grid point to the next, one jumps between these two solutions, and that 
is equivalent to a high frequency mode. 



85 



B Numerics 



B.3.2 Elliptic solvers 

Solving elliptic equations is expensive. The elliptic solver needs the majority (more than 90%) of 
the total computing time in the Tiger code. This is partly my fault, as 1 use methods that are known 
to be subideal and very s low. 

Currently 1 use PETSc jBBG+0lllBGMS01llBGMS97l to solve the elliptic equations. PETSc uses 
a Newton-like method to reduce the nonlinear equations to linear ones, and then uses Krylov- 
subspace methods to solve these. 1 calculate the Jacobian that is necessary for the Newton-like 
method numerically. This is known to be a very expensive operation; it is about ten to twenty 
times slower than a hand-coded Jacobian. On the other hand, hand-coded Jacobians would need 
to be written and debugged. 

1 also have a simple iterative Jacobi solver for testing purposes. This solver is only a toy, and is 
only usable for very small problems. 

1 estimate that a nonlinear multigrid solver would be several orders of magnitude faster than 
the way 1 (ab)use PETSc, as it does not need to calculate the Jacobian. An experimental implemen- 
tation of a nonlinear full approximation storage multigrid solver for the two-dimensional Poisson 
equation with excision shows promising results. 



B.4 Coding equations 

As mentioned in section 17.11 1 chose to code the equations by hand instead of using a symbolic 
algebra package. This only makes sense because it is possible to make the computer code look 
sufficiently close to an explicit mathematical notation. The standard mathematical notation used 
in general relativity is full of context-sensitive abbreviations, such as the Einstein sum convention, 
that cannot be used in an unambiguous notation. 

One of the key points to obtain a readable source code is to make use of the local character of the 
equations, leading to a pointivise programming style. A tensor such as Kjj is really a tensor field 
Kij{x^), which becomes a grid function Kij{x,y,z) after discretisation. Dealing with a set of five 
indices (z, x,y,z) at the same time leads to an awkward notation. Hence 1 chose to deal with one 
grid point at a time only. 

As an example, 1 present the code for one of the more complicated equations, namely the time 
evolution equation of the extrinsic curvature. It is given in appendix lA. 1 .2l as: 



Rjj — IKjiK^j + KKjj 
K,jl5l, + Kal3lj + l3'K,j^i 



(B.20) 



When translated into a Fortran subroutine for Cactus, this equation becomes 



subroutine calc_extcurv_dot k 

(gu, ri , kk, dkk, alfa, ggalfa, beta, dbeta, kk_dot) 



CCTK_REAL, 
CCTK_REAL, 
CCTK_REAL, 
CCTK_REAL, 
CCTK_REAL, 
CCTK_REAL, 



intent (in) 
intent (in) 
intent (in) 
intent (in) 
intent (in) 
intent (out) 



gu(3,3) 
ri(3,3) 

kk(3,3), dkk(3,3,3) 
alfa, ggalfa(3,3) 
beta(3), dbeta(3,3) 
kk_dot(3,3) 



86 



B.4 Coding equations 

integer : : i, j ,k,l 

! K_ij,t = - D_i D_j alpha 

! + alpha [ R_ij - 2 K_ik K~kj + K K_ij ] 

! + K_kj d_i beta~k + K_ik d_j beta~k + beta~k d_k K_ij 

do i=l,3 
do j=l,3 

kk_dot(i,j) = - ggalfa(i,j) + alfa * ri(i,j) 
do k=l,3 
do 1=1,3 

kk_dot(i,j) = kk_dot(i,j) & 

- 2*alfa * gu(k,l) * kk(i,k) * kk(j,l) & 
+ alfa * gu(k,l) * kk(k,l) * kk(i,j) 

end do 

kk_dot(i,j) = kk_dot(i,j) & 

+ kk(k,j) * dbeta(k,i) + kk(i,k) * dbeta(k,j) & 
+ beta(k) * dkk(i, j ,k) 

end do 
end do 
end do 

end subroutine calc_extcurv_dot 

In the above routine, the variables gu, ri, etc. represent the corresponding tensors at a single 
point in the time slice only. This routine is called once for every grid point for every right hand 
side evaluation. Such a coding convention needs powerful compilers and optimisers to yield good 
performance, and likely cannot be used on vector machines. This also shows that programming 
styles are heavily influenced by the current software and hardware technology. However, I strove 
for clarity and simplicity rather than for performance. 
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C Definitions 



C.l Glossary of terms 

Agave: A numerical code for the vacuum Einstein equations, descending from the Grand Chal- 
lenge 

Cactus: A framework f or solving three-dim ensional time-dependent systems of partial differen- 
tial equations, see lABD+OlllABG+OlllCacl 

coordinate condition: A condition that determines lapse a and shift j3' (equivalent to the term 
gauge condition) 

excision: A method to treat singularities. One cuts a region of the spacetime that contains the 
singularity out of the simulation domain, introducing an inner boundary. This is possible 
as long as the boundary is located within the event horizon. As no physical information 
can propagate out of the event horizon, the details of how the excision boundary is handled 
cannot influence the simulation outside the event horizon 

gauge condition: A condition that determines lapse a and shift /3' (equivalent to the term coordi- 
nate condition) 

gauge evolution condition: A choice of lapse a and shift /3' that imposes no restrictions on the 
three-metric y,y and extrinsic curvature Kij. For example, geodesic slicing (a = 1) is a gauge 
evolution condition 

gauge fixing condition: A condition on the three-metric y,y and extrinsic curvature Kij that also 
leads to a prescription for lapse cc and shift /3' via the ADM time evolution equations. For 
example, maximal slicing (K = 0) is a gauge fixing condition 

IVIaya: A numerical code, descending from Agave, see IPenI 

performance: (of a code) the quality of the result of a simulation, i.e. a measure of the magnitude 
of the numerical errors. This has nothing to do with the speed of the code on a particular 
hardware 

punctures: A method to treat singularities. The singularity is contained in a conformal factor t/j, 
and all other quantities can be chosen so that they remain finite. This is the case e.g. for Brill- 
Lindquist black holes (see section l6.3.H . The singularity is placed in between grid points. The 
conformal factor and its derivatives are known analytically. With maximal slicing {K = 0) 
and normal coordinates (/3' = 0), the conformal factor is then constant in time 

robust stability: A practical definition for the notion of stability of a code that can be tested em- 
pirically See ISGBW00..SSW02J 
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C Definitions 



Tiger: A numerical code, descending from Maya 



time slice: Another term for a spacelike hypersurface 



C.2 Abbreviations 

ADM: Arnowitt-Deser-Misner. An initial-boundary-value formulation for the Einstein equa- 



tions, see I ADM62J 
AEI: the Albert-Einstein-Institut in Golm near Potsdam, Germany 
BH: black hole 

BSSN: Baumgarte-Shapiro-Shibata-Nakamura. An ADM-like initial-boundary-value formula- 
tion for the Einstein equations, see LNOK87..N089..SN95llBS99l 

ctADM: conformal traceless ADM system 

KST: the Kidder-Scheel-Teukolsky system; see iKST^OOl 

MPI: Max-Planck-Institut, e.g. the AEI, which is also called the Max-Planck-Institut fiir Gravi- 
tationsphysik 

MTW: the book Gravitation by Misner, Thorne, and Wheeler; see IMTW73I 

NS: neutron star 

PSU: Perm State University 

TAT: Theoretische Astrophysik Tubingen, the Abteilung fiir Theoretische Astrophysik in the In- 
stitut fiir Astronomie und Astrophysik of the Fakultat fiir Mathematik und Physik at the 
Eberhard-Karls-Universitat Tiibingen (yes, this is in Germany) 

TGR: the name of the Tiibingen General Relativity (Tiger) code 

C.3 Symbols 

It goes without saying that Greek letters denote four-vector indices ... 3, while Latin letters de- 
note three-vector indices 1 ... 3. I follow the sign conventions of MTW IMTW73I . 

Four-quantities: 

rin-v Minkowski metric, rj^y = diag( — 1, +1, +1, +1) 

g^y four-metric 

Gfxy Einstein tensor 

Tf^y energy density tensor 

Three-quantities: 
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C.3 Symbols 



5ij identity tensor. <5,y = diag(+l, +1, +1) 

Yij three-metric, see appendix lA.l.ll 

Kjj extrinsic curvature, see appendix lA.1.21 

a lapse, see appendix lA.l.ll 

^' shift, see appendix lA.l.ll 

Fyj. connection coefficients, see appendix lA.1.31 

r' contracted connection coefficients. V = T^'^Tj^. 

Rjj Ricci tensor, see appendix lA.1.31 

H Hamiltonian constraint, see appendix lA.1.41 

Mi momentum constraint, see appendix lA.1.41 

Conformal quantities: 

~ a tilde denotes a conformal quantity, which has its indices raised and lowered with the 
conformal metric 

xp i/)^ is the conformal factor, see appendix lA.2.11 

Yij conformal three-metric, see appendix lA.2.11 

K trace of the extrinsic curvature, see appendix IA.2. II 

Aij traceless extrinsic curvature, see appendix IA.2. II 

Aij conformal traceless extrinsic curvature, see appendix IA.2. II 

Vi vector potential for the extrinsic curvature, see appendix lA.2.41 

Coordinate quantities: 

" a bar denotes a coordinate quantity, which has its indices raised and lowered with the 
coordinate metric (i.e. for Cartesian coordinates) 

h-ij traceless conformal three-metric, see appendix lA.2.51 

F, gauge condition for the metric, akin to a constraint variable, see appendix lA.2.51 
W, vector potential for the traceless three metric, see appendix lA.2.61 
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And crawling on the planet's face 
Some insects called the human race 
Lost in time, and lost in space 
And meaning. 



